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COMPUTER  MODEL  PREDICTING  THE  THERMAL  RESPONSE 


5 TO  MICROWAVE  RADIATION  OF  A SIMULATED  BIOLOGICAL  STRUCTURE 


r 


\|^'  1.  INTRODUCTION 

■'ills  paper  describes  a method  of  computing  the  thermal  response  of  an 
autothermal 1y  regulated  body  such  as  a biological  body  to  a source  of  micro- 


wave radiation.  The  description  of  this  "wetiTod  is  divided  into  five  parts. 
It  includes  (1)  a discussion  of  the  symbols  (and  their  units)  used  in  develop- 
ing the  solution,  (2)  the  induced  electromagnetic  field  distribution  and  the 
power  density  distribution  that  represents  the  source  term  for  the  heat  trans- 
fer problem,  (3)  the  modified  heat  operator  eigenvalues  and  eigenfunctions 
associated  with  a Newton  cooling  law  boundary  condition,  (4)  the  computation 
of  temperature  excursions  induced  by  microwave  radiation  including  complex 
pulse  heating  schemes,  and  (5)  a discussion  of  the  spreading  of  temperature 
distributions  with  time  in  three  types  of  simulated  biological  structures. 
These  are  discussed  in  Sections  3. 1-3. 5 respectively. 

In  [2]  we  developed  a computer  model  to  determine  the  temperature  distri- 
bution in  a penetrable,  homogeneous,  and  spherically  symmetric  body  that  has 
been  irradiated  by  microwave  radiation.  Heat  removal  by  blood  flow  could  be 
considered,  provided  that  only  one  layer  was  used  in  the  model.  In  the 
present  paper  a shooting  method  for  solving  the  eigenvalue  and  eigenfunction 
determination  problem  for  a multilayered,  penetrable,  but  spherically 
symmetric  scatterer  is  solved  when  the  heat  equation  describing  the  microwave 
heating  includes  the  possibility  of  blood-flow-heat-removal  terms  in  some,  but 
not  necessarily  all  layers.  This  innovation  is  described  in  Section  3. 


OJ  - 
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Also,  originally  our  program  in  [2]  experienced  some  difficulty  in  computing 
expansion  coefficients  used  in  determining  the  induced  electric  field  when  the 
frequency  of  the  incoming  radiation  was  high  (>10  GHz)  or  when  the  radius  of 
the  outer  sphere  was  as  large  as  48  cm;  the  procedure  by  which  we  overcame 
this  difficulty  is  described  in  Section  2.  Some  experimental  microwave 
bioenvironmentalists  look  for  a nonthermal  microwave  effect  and  consequently 
attempt  to  control  temperature  in  their  microwave  exposure  systems  by  using  a 
complex  microwave  pulse  heating  scheme  with  a low  duty  factor.  Section  4 
therefore  contains  a description  of  a formula  which  permits  one  to  express  the 
expansion  coefficients  associated  with  a complex  temporal  heating  pattern  in 
terms  of  integrals  with  respect  to  only  the  spatial  variables.  Finally  we 
note  that  several  people— including  MacLatchy  and  Clements  [9]  and  Washisu  and 
Fukai  [11]  have  proposed  microwave-induced  temperature  excursions  as  a 
nonperturbing  method  of  measuring  or  estimating  field  strengths. 
Consequently,  in  Section  3.5  of  this  paper  we  have  included  a discussion  of 
the  potential  and  limitations  of  this  method  of  field  measurement.  Also, 
because  microwave  heating  may  be  used  to  treat  tumors  in  humans  (c.f.  Zimmer 
et  al.  [12]),  we  give  in  Section  3.5  some  new  computer  calculations  showing 
possible  thenal  effects  on  simulated  biological  structures. 
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2.  RESULTS  AND  DISCUSSION 


In  this  paper  the  authors  extend  the  computer  model  which  generated  the 
results  of  [2],  We  consider  as  before  that  a plane  wave  irradiates  a spheri- 
cally symmetric  structure  in  the  manner  indicated  in  Figure  2,1.  We  allow, 
however,  time  profiles  similar  to  that  of  the  PAVE  PAWS  radar  so  that  heating 
from  any  radar  emission  can  be  estimated  directly.  We  note  that  with  this 
capability  and  the  possibility  of  estimating  temperature  derivatives  that  one 
can,  by  solving  the  equations  of  thermoelasticity,  describe  radar  acoustic 
effects  in  a quantitative  way.  This  is  important  in  view  of  large  efforts  by 
other  branches  of  the  armed  services  to  study  this  effect.  We  can  also,  by 
going  to  more  general  geometries  and  using  an  integral  equation  method, 
describe  quantitatively  the  effect  of  microwave  radiation  on  biochemical  pro- 
cesses and  fetal  development  which  are  strictly  thermal  in  nature. 

Figures  2.2-2.10  compare  computer  model  predictions  with  measurements 
made  by  John  G.  Burr  at  Brooks  AFB  and  give  a comparison  of  our  ability  to 
predict  spatial  variations  in  temperature  for  two  radiofrequencies  (RF),  1.2 
GHz  and  2.5  GHz,  and  for  a short  30-s  and  for  a longer  3-min  exposure.  The 
capability  of  predicting  the  thermal  response  to  pulsed  radiation  is  demon- 
strated in  Figure  2,11. 

We  now  discuss  the  effect  of  blood  flow  in  removing  heat  from  a living 
system  subjected  to  microwave-radiation-induced  thermal  excursions.  John  Burr 
and  Jerome  Krupp  of  the  Radiation  Sciences  Division  of  the  USAF  School  of 
Aerospace  Medicine  reported  in  [3]  the  results  of  Figure  2,12  showing  tempera- 
ture measurements  in  the  head  of  a living  and  dead  Macaca  mulatta.  In  Figure 
2.13  we  use  blood  flow  rates  supplied  in  [9]  to  estimate  the  effect  of  the 
blood  flow  term  in  giving  a lower  predicted  value  of  a radiation-induced  tem- 
perature increase. 
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Figure  2.2.  Temperature  rise  along  the  z-axis  of  a 3.3-cm  radius  homogeneous 
muscle-equivalent  sphere  exposed  to  1.2  GHz,  continuous  wave 
(CW),  70  mW/cm^,  RF  in  the  far  field  for  30  s. 
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Figure  2,3.  Temperature  rise  along  the  x-axis  of  a 3,3-cm  radius  homogeneous 
muscle-equivalent  sphere  exposed  to  1.2  GHz,  CW,  70  mW/cm^,  RF  in 
the  far  field  for  30  s. 
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Figure  2.4.  Temperature  rise  along  the  y-axis  of  a 3.3-cm  radius  homogeneous 
muscle-equivalent  sphere  exposed  to  1.2  GHz,  CW,  70  mW/cm^,  RF  ir, 
the  far  field  for  30  s. 
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Figure  2.5.  Temperature  rise  along  the  z-axis  of  a 3.3-cm  radius  homogeneous 
muscle-equivalent  sphere  exposed  to  2.5  GHz,  CW,  100  mW/cm^,  RF 
in  the  far  field  for  30  s. 


-3  -2  -I  0 +1  +2  +3 

Y -AXIS  (cm) 


Figure  2. 


. Temperature  rise  along  the  y-axis  of  a 
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Figure  2.8.  Temperature  rise  along  the  z-axis  of  a 3.3-cm  radius  homogeneous 
muscle-equivalent  sphere  exposed  to  1.2  GHz,  CW,  70  mW/cm^,  rf  in 
the  far  field  for  3 min. 


Figure  2. 


. Temperature  rise  along  the  x-axis  of  a 3.3-cm  radius  homogeneous 
muscle-equivalent  sphere  exposed  to  1.2  GHz,  CW,  70  mW/cm^,  RF  in 


the  far  field  for  3 min 


Figure  2.10,  Temperature  rise  along  the  y-axis  of  a 3.3-cm  radius  homogeneous 
muscle-equivalent  sphere  exposed  to  1.2  GHz,  CW,  70  mW/cm^,  RF  in 
the  far  field  for  3 min. 
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2.11.  The  predicted  and  measured  temperature  excursion  versus  time  at 
the  center  of  a 3.3-cm  radius  homogeneous  muscle-equivalent 
sphere  at  1.2  GHz,  CW,  70  mW/cm^. 
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Figure  2.12.  Temperature  excursion  in  the  midbrain  of  a living  (blood  flow 
case)  and  dead  (no  blood  flow  case)  Macaca  mulatta  (rhesus 
monkey)  head  exposed  to  70  mW/cm^,  CW,  RFR  in  the  far  field,  1.2 


We  see  that  there  are  qualitative  differences  in  the  curvature  of  the 
temperature  versus  time  curves  in  Figures  2.12  and  2.13  that  are  at  this  point 
unexplained  although  the  measured  and  predicted  values  seem  to  be  reasonably 
close. 

Finally  we  give  a comparison  between  our  computer  model  and  the  simpler 
model  developed  by  Kritikos  and  Schwan  [5].  This  comparison  is  given  in 
Figures  2.14  and  2.15. 

We  note  that  consideration  of  the  thought  experiment  depicted  in  Figure 
2.16  makes  it  obvious  that  there  is  such  a thing  as  a nonthermal  effect.  We 
consider  a structure  that  will  not  drift  in  a microwave  field  but  which  will 
necessarily  respond  differently  to  a microwave  field  and  to  an  equivalent 
amount  of  thermal  energy.  We  consider  a simple  molecule  with  three  charges  in 
a row  connected  by  two  chemical  bonds  that  we  approximate  by  two  linear 
springs  with  identical  spring  constants.  The  outside  two  moieties  have  charge 
q and  the  middle  moiety  has  charge  -2q, 

Since  all  three  masses  are  the  same,  the  thermal  energy  of  the  solvent 
will  act  in  the  same  way  on  all  three  moieties,  but  the  electric  field  will 
exert  twice  as  much  force  on  the  inner  moiety. 


COMPRRISON  OF  KRITlKOS-SCHWfiN  AND 
MIE  SOLUTION  OENERHTED  SOURCE  TERM 


Figure  2,14.  Comparison  of  the  Kritikos  and  Schwan  source  term  used  in  [5]  and 
the  Mie  solution  generated  source  term  used  in  this  paper.  The 
magnitude  of  the  Kritikos  and  Schwan  source  term  is  10,000  W/m^. 
The  Mie  solution  assumes  that  the  incident  power  is  10  mW/cm^ 
(field  strength  * 194,09  V/m),  that  the  frequency  is  1000  MHz, 
that  the  real  part  of  the  relative  permittivity  is  34.4,  that  the 
ionic  plus  polarization  current  conductivity  = o'  + ueoe"  = 0,8 
mhos/m,  and  that  the  outer  boundary  of  this  scattering  body  is  a 
sphere  whose  radius  is  5 cm. 
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Figure  2.15.  Comparison  of  the  Kritikos-Schwan  predictions  in  [5]  (marked  with 
an  *)  and  our  solution  (smooth  curve).  We  assumed,  following 
Kritikos  and  Schwan,  that  the  blood  flw  was  normal  (b  = 0.00186 
cal/cm^/s)  and  that  the  exposure  time  was  200  s;  we  used  the 
parameters  K = 0.001  ca1/cm/®C,  p = 1.0  g/cm^,  and  c = 1.0 
cal/g  • °C  that  were  used  in  [5]. 


> FIELD  LINES  PAHALLEL 
TO  LIWE  JOININfr  THE 
CENTERS 


Figure  2.16.  Electromagnetic  field  interaction  model  for  which  there  would  be 
a nonthermal  effect. 
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3.  MATHEMATICAL  PRELIMINARIES 


3.1.  Notation 

The  variables  used  in  this  paper  are 


ENGLISH 

A^  = dimensionless,  coefficient  of  the  radial  eigenfunction 
that  is  regular  at  the  origin. 


expansion  coefficient  for  the  odd  regular  vector  wave 
functions. 


a.('’’’”){t)  = the  temperature  decay  factor  of  the  solution  u associated 
k 

with  the  radial  eigenvalue  and  the  Legendre 

transform  lH!  > 
n 


b(r)  » the  product  of  the  number  of  grams  of  blood  per  gram  of 
tissue  per  second,  the  tissue  density  in  grams  of  tissue 
per  cubic  centimeter  of  tissue,  and  the  specific  heat  of 
the  blood  (typically  b = .0122), 

= dimensionless  coefficient  of  the  radial  eigenfunction 
that  is  singular  at  the  origin. 


b(^  p)  ' expansion  coefficient  for  the  even  regular  vector  wave 
functions, 

b(m,n)(t)  3 temperature  decay  factor  of  the  source  term  S 

(associated  with  the  radial  eigenvalue  and  the 

Legendre  transform  , 


c(r)  = tissue  specific  heat  in  calories  per  gram  degree 
centigrade  (typically  c = .84), 

t = the  electric  field  intensity  in  volts  per  nieter  (10 

milliwatts  per  square  centimeter  corresponds  to  194.087 
volts  per  meter), 

f = frequency  in  Hertz, 

H = Newton  cooling  constant  in  calories  per  square  centimeter 
per  degree  per  second  (typically  H ^ .0000572), 

fi  = magnetic  field  intensity  in  Henrys  per  meter  (10  milli- 
watts per  square  centimeter  corresponds  to  .5151  Henrys 
per  meter) 

h^  * spherical  Hankel  function  “ 0^  - iYp  that  is  used  in 
expanding  the  electromagnetic  fields  in  Tesseral 
harmonics, 

^(^(X,r),r,n)  * the  radial  eigenfunction,  used  in  expanding  the 
temperature,  that  is  nonsingular  at  the  origin, 

= the  spherical  Bessel  function  of  order  n, 

K = thermal  conductivity  in  calories  per  centimeter  per 
degree  centigrade  per  second  (typically  K = .0012), 

m = the  index  of  the  finite  cosine  transform  (m  = 0 or  1 in 
our  application). 
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n = the  index  of  the  Legendre  polynomial  used  in  expanding 
the  field, 

P^(cos(0))  = the  associated  Legendre  polynomial, 

r = the  distance  from  the  center  of  the  scatterer  in  centi- 
meters, 

= the  rudius  of  the  ith  bounding  sphere  in  centimeters, 

S = the  source  term  for  the  heat  equation  in  calories  per 
cubic  centimeter  per  second, 

S(X,r)  = (Xp(r)c(r)  - b(r))/K(r). 

^^(x)  = a constant  value  of  ^{X,r)  occurring  when  < r < R^ , 
t = time  in  seconds, 

u = temperature  excursion  above  the  ambient  temperature, 
y^  = spherical  Bessel  function  of  the  second  kind, 

^n+1/2  " order  Bessel  function  of  the  second  kind  (Weber 

function  of  order  n+1/2), 

X(l(x,r),r,n)  = the  radial  eigenfunction,  used  in  expanding  the 
temperature,  that  is  singular  at  the  origin, 

Z(^,k)(r‘)  = the  radial  eigenfunction  associated  with  the  eigenvalue 


GREEK 


“(^P) 


= expansion  coefficient  for  the  odd  singular  vector  wave 
functions, 


{X,R,ri)  = (X),r,n), 

a^(X,R,n)  = K^[(3/Sr)^(^.(x),r,n)]  evaluated  at  r = R, 


p)  = expansion  coefficient  for  the  even  singular  vector  wave 
functions. 


6.(X,R,n)  = Y(S.(X),R.n], 

i^(X,R,n)  = J (a/3r)V^(^^(X),r,n)]  evaluated  at  r = R, 
^^•(x,R^,n)  = + »^)  " »f')» 


e = permittivity  in  farads  per  meter. 


0 = spherical  coordinate— angle  of  ray  to  a point  with  the 
positive  z-axis. 


X = eigenvalue  associated  with  radial  harmonics, 
p = density  in  grams  per  cubic  centimeter. 


a = conductivity  in  ohms  per  meter, 

T = dummy  variable  of  integration  used  in  expressing 

temperature  decay  factors  as  a convolution  integral, 

<j)  = spherical  coordinate  of  the  x-y  plane. 


0)  = frequency  in  radians  per  second. 
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x(Np,T,Tp,Tp)  = a cutoff  function  for  the  temporal  envelope  of  the  pulse 
heating  scheme. 

MISCELLANEOUS 
3 = partial  derivative  symbol 
ENGLISH  SCRIPT 

pulse  heating  scheme  temporal  envelope  function, 

C * the  finite  cosine  transform, 

m * 

= the  Legendre  transform, 

^(Tp»Tp.Np.T^,T)  « part  of  a temperature  decay  factor  associated  with  a 

heating  pattern  defined  by  equation  (3.4.14), 

Jp»NpJ^j»T)  = part  of  a temperature  decay  factor  associated  with  a 
complex  pulse  heating  scheme  defined  by  equation 
(3.4.15), 

and 

j.(n  ® the  radial  transform  used  in  getting  expansion 

coefficients  to  express  a function  of  r in  terms  of 
radial  eigenfunctions  that  satisfy  the  Newton  cooling  law 
boundary  condition. 

Other  notation  that  is  introduced  in  the  text  of  the  paper  is  defined  and  used 
locally. 


3.2.  Induced  Electromagnetic  Field  Distribution 


A new  practical  method  of  developing  Tesseral  harmonic  expansion 
coefficients  for  the  electromagnetic  field  induced  in  a penetrable  scatterer 
with  spherical  symmetry  is  described  here.  Our  numerical  technique  will 
permit  us  to  use  the  Mi  e-solution  method  to  determine  the  response  of  a body 
with  a large  size  to  a higher  frequency  radiation  than  we  could  with  the 
standard  methods  described  in  the  references  of  [1]. 

The  electric  field  induced  in  the  pth  interior  region  by  a wave  of  the 

form 


= E^exp  (-iu) 


oCt- 


(3.2. 


is  given  in  the  pth  region  by 


ill  (^»P)  (1,4) 


ib  n"-‘>  . a 

(£,P)  (l,Jt)  (t,p)  (l,t) 


- iB 


•*-(e,3) 

N 

(A.P) 


(3.2. 


where 


M 


(ej) 


(i.n) 


^^si(kpr)p‘(cos(e))stnU)?, 


- Zn(kpr)(d/d0)(Pn(cos(e)))cos(*)e^,  (3.2, 


5(o.J) 

(i.n) 


■ Zn(kpr)(d/d0)P^{cos(e))sin(4))e^, 


{3.2.4) 


and 


where 


"Ir.n)  = 


+ (^) (9/3") r) ) (d/de ) (P^ (cos (6 ) ) )cos (^)eQ 

kpP 

■ ((k  r)lin(6)) 


Sjo.j)  = 2J(k  r)p‘(cos(8))sinU)? 

(i,n)  kpP  n p n r 

+ ~ (3/9r)(rz^(kpr))(d/d9)(pJ(cos(0)))sin(^)eg 

^ (;k'~r')in'C?TJ  (3/3>-)(rzj(kpr))pJ(cos(e))cos(*)l^.  (3.2.6) 


kp  ^ sign  (u) 


.2  + jjjj~rjj3- 


n 


-yew 


2 + /vW  * 


]■ 


(3.2.7) 
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In  (3.2.3)  - (3.2,6)  functions  are  the  associated  Legendre  polynomials 
and 


j = 3 

zj(z)  = , (3.2. 

j^{z)  if  j = 0 

v/here  and  are  respectively  the  spherical  Bessel  functions  of  the 
first  and  second  kind.  Part  of  the  difficulty  is  that  we  cannot  use  (3.2.7) 
to  compute  h^  even  if  we  know  and  y^  exactly.  For  example,  z = u + iv 

implies  that 

hj(z)  = (l/z)Csin(u)(cosh(v)  - sinh(v)) 

+ i cos(u)(sinh(v)  - cosh(v))]  (3.2. 

which  is  uncomputable  on  a digital  computer  if  v is  large  enough  so  that 
cosh(v)  and  sinh(v)  are  indistinguishable, 

A better  way  is  the  use  of  the  formula 


hj(z)  = i""'^z"^exp(iz)  I (n+l/2,k)(-2iz)'^, 

k=0 


(3.2.1 


coupled  with  the  Hankel  symbol  formula. 


(n+l/2,k) 


(n+k)!  _ (n+k)(n+k-l)'»»» (n+l)n(n-l)»»» 
ki(n-k)!  k! 

k ,(n-(k-2i))(n-(k-(2i-l))) 
n ( }, 


(n-(k-l)) 

(3.2.1 


when  the  complex  number  z is  such  that  (n+l/2,k)(-2iz)  are  of  such  a 
size  that  round-off  error  is  not  encountered  in  the  computation  of  (3.2.10). 
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The  reader  can  verify  (3.2.10)  by  induction  using  the  three-term  recursion 
formula 


h‘ . = h'  - h‘ 

n+1  2 n n-i 


(3.2.12) 


is  satisfied  and  by  showing  that  (3.2.10)  is  true  for  n = 0 and  n = 1 


since  from  (3.2.8)  we  know  that 


h‘(2) 


(3.2.13) 


and 


hj(z)  = (• 


sin(z)  cos(z).  cos(z)  sin(z) 
- 1 + il-n ■ ■ ■ - 


-).  (3.2.14) 


For  intermediate  values  of  z,  another  method  must  be  sed  to  compute 

1 

h^(z).  We  observe  that  equation  (3.2.12)  implies  that 


1 


1 


(3.2.15) 


The  basic  idea  is  to  write 


n+1/2 


= = ^(n-i-H/2)^^^ 

^n^^^  ^(n+1/2)^^^ 


(3.2.16) 


and  then  observe  that  equations  (3.2.15)  and  (3.2.16)  imply  that 


'(n+1/2) 


2(n+l)-l 


'(n+1+1/2) 


(3.2.17) 
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Hence,  v = n + 1/2  and  n = v - 1/2  and  equation  (3.2,17)  imply  that 


5 


ri 


N 


a 


V 


1 

3v+l 


(3.2.18) 


Thus,  from  (3.2,16)  we  get  immediately  a continued  fraction  expansion 


a 

V 


2v 

z 


1 

2(v-H)  . 1 
" ^+2 


(3,2.19) 


et  cetera,  which  by  the  following  Lemma  is  always  convergent. 

Lemma  (Wall  [10],  p.  50).  We  have  uniform  convergence  of  the  continued 
fraction. 


c « 


”1  * 


“■a* 


b + !j 

3 b + 


if  there  exists  constants  gpe(0,l)  such  that 


p+i 


b b 
P P+i 


i - 9p)Vl' 


In  our  situation 


b 


0 


z 


b =* 

P z 


P “ l-»2,... 


(3.2.20) 


(3.2.21) 


(3.2.22) 


(3.2.23) 
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and 


3 


N 


Si 


I*. 


(3.2.24) 


Thus,  for  every  z there  is  an  > 0 such  that  if  p > then 


P+i 


b b 
P P+i 


4(v+p)(v+p+l) 


.1  ^^"9p)Sp+i  (3.2.25) 


provided  that  is  such  that  (N^+v)  jzj  and  Qp  = 1/2  for  all  p.  Thus 
in  view  of  the  Lemma  the  continued  fraction  expansion  theoretically  converges 
for  all  z ♦ 0. 

The  idea  then  is  not  to  compute  the  spherical  Bessel  functions  of  the 
second  kind  y^  at  aP,  but  rather  use  a direct  method  for  obtaining  the  hj^. 
Observe  that 


h*  = 


% K , h (-i)exp(iz) 

^ i_  ••  * * 


(3.2.26) 


h , h _ 
n-i  n-2 


The  functions  J^(z)  used  in  the  expansion  are  successfully  computed  by  the 
methods  of  Lentz  [8]. 
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3.3.  Heat  Operator  Eigenvalues  and  Eigenfunctions 
for  a Nevfton  Cooling  Law  Boundary  Condition 


3.3.1.  The  Radiative  Heat  Transfer  Problem.  From  the  E field  determina- 
tion of  the  preceding  section,  we  develop  an  expression  for  a source 

S - div(E  X H)/(10^  X 4.184)  (3.3.1) 

of  internal  energy  generation  which  is  used  as  a term  in  the  heat  equation, 

pc  — - div(K  grad(u))  + bu  = S,  (3.3.2) 

3t 

where  pc  is  the  product  of  density  and  specific  heat,  b is  a blood-flow 
cooling  term,  and  K is  the  thermal  conductivity.  Assume  that  the  scattering 
body  is  a union  of  material  regions  bounded  by  spheres  r = for  i in 
{1,...,N}  (with  N _<  6 in  our  computer  program)  and 

0 * Rq  < Ri  < ...  < R|^,  (3.3.3) 

Assume  that  p(r),  c(r),  and  K(r)  have  the  constant  values  p. , c^, 
respectively  for  R^._^  < r < R^.  Then  for  R^_^  < r < R^.  equation  (3.2)  may 
be  written 


p^c^(a/3t)u  = 


3r 


3r 


r sin(9)  39  39 


+ 


1 3^Ui 

— — J 

sin(9)  3(|i' 


- b^u  + S, 


(3.3.4) 


p 

8 


i 


I 


where  the  initial  condition  is  that 

u(r,0,*,O)  = 0, 


(3.3.5) 


continuity  of  temperature  and  heat  flux  implies 


Lim  u(Rj-e,0,({>,t)  = Lim  u (Rj+e,0,(|»,t)j 
e 0 ^ e ^ 0 ^ 


(3.3.6) 


and 


Lim  K. (3/3r)u(R.-e,0,({i,t)  - Lim  K.  (3/3r)u(R.+e,0,(|),t),  (3.3.7) 

e 0 ^ ^ e - 0 ’ 


and  the  Newton  cooling  law  implies  that 

K^(3/3r)u(Rf^,0,^,t)  + Hu(Rj^,0,t,t)  « 0. 


(3.3.8) 


We  define  the  finite  cosine  transform  of  the  temperature  excursion  u(r,0,(ii,t) 
by  the  rule, 


(<^m‘^){f»6,t)  = (I/tt)/  u(r,0,^,t)cos(m<fr)d(fr, 

I" 


(3.3.9) 


for  positive  integers  m and 


(<7  u)(r,0,t)  = (l/2it)/  u(r,0,^,t)d4.. 

^ -TT 


(3.3.10) 


We  define  the  Legendre  transform  operator  on  the  temperature  excursion 

If 


u(r,a,(j,t)  by  the  rules, 


(i'nU)(r,<!),t)  = 


(^)(-|^^^]A(r,e,^,t)Pj’(cos(0))sin(0)d0 
1 (n+m)*  0 '' 


(3.3.11) 
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= (il!ll)/^u(r.6,(fr,t)P  (cos{0))s1n(0)d9. 
n 10  " 


(3.3.12) 


Thus,  If  we  combine  (3.3.9),  (3.3.10),  (3.3.11),  (3.3^12),  and  (3.3.4) 
that  if  p,  c,  and  K are  simply  functions  of  r,  then 

(3/3t)i;;e^u  = (J_)((3/3r)(r^K(r))(3/3r)I^c^u 
per 

- K(r)r.{ntl)£Uc^ul  - (b/(p))£™C^u  + £|;(7^(S/pc) 


Let  us  attempt  to  write 


and 


where 


k»l 


£:(:„(S/PC)(r.t)  = i b<"-">(t)Z,  ,)(r). 
k=l 

(d/dr)(K(r)r^(d/dr))Zjj,^l^j(r)  + 


Limy 

e-»-0^(n,k) 


Limy 

e-^^(n,k) 


we  see 

(3.3.13) 

(3.3.14) 

(3.3.15) 

(3.3.16) 

(3.3.17) 

(3.3.18) 


for 


where  the  k)  ^re  positive  numbers  for  which 

Thus,  we  conclude  that 

and  consequently  that 


By  defining  for  every  function  g(r) 


■(n.k) 


(g) 


^0 


(3.3.19) 


(3.3.20) 


(3.3.21) 


(3.3.22) 


we  see  that 


3.3.2.  Eigenvalue  Determination.  We  wish  to  describe  here  a computer 
algorithm  for  computing  the  and  the  k)^^^  ^ 

are  nonnegative  piecewise  constant  functions.  While  this  is  trivial  when  b(r) 
is  a constant  function,  the  problem  becomes  interesting  when  b(r)  is  posi- 
tive  in  some  intervals  and  is  identically  zero  in  others. 
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To  begin  with  we  define 


r)  - ^p(r)c(r)  - b(r) 
~ ’ K(r) 


(3.3.24) 


and  we  obtain  in  each  interval  op®  three  different  classes  of 

solutions  of  the  problem  (3.3.16)  - (3,3.19)  depending  on  whether  or  not 
^(X(n  is  uniformly  positive,  zero,  or  negative  in  this  interval.  We, 

therefore,  write 

Z(r,X)  = A^jJ(^(X,r),r,n)  + B^Jf(^(X,r),r,n)  (3.3.25) 

and  require  that 


L,„  AjJ{S(i,r).r,n) 
r*0 r ■ ‘ 


(3.3.26) 


Bi=0 


(3.3.27) 


In  general  for  ^(X,r)  > 0 we  have  setting 


2 = r /h^  = r/l(^ 


(3.3.28) 


the  fact  that  (3.3.16)  is  equivalent  to 


r^(d/dr)^Z^  + 2r(d/dr)Z^  + r^S{X,r)Zj^  - n(n+l)Z^  = 0 (3.3.29) 


and  if  we  change  variables  by  the  rule 


z = /S(X,r)r 


(3.3.30) 
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and  observe  that 


we  see  that  if 


d ^ dz  d 
dr  dr  dz 


then 


z^W"  + 2zW'  + (z^-(n(n+l)W  = 0, 
which  implies  that  in  R^)  we  have 


W{z)  = Aijn(z)  + Biyn(z), 


where 


(3.: 


{3.: 


(3.: 


(3.; 


j„u) 


n (2k-l) 

k=l 


maO  ?)  \ 

ml  n (2k4l+2n) 
•"•kaO 


(3.: 


Thus,  to  make  (3.3,34)  consistent  with  (3.3.26)  when  i = 1,  we  must  have 

n+1  

A = l/((  It  (2k-l))(/S(X.r))''|.  (3. 

‘ k’l 

If  ^(A,r)  = 0 we  have 


J(S(X,r),r,n)  = r” 


(3. 


and 


n 


K'-'- 

.*  N 


I(S{X,r),r,n)  = r 


-n-i 


(3.3.38) 


If  S(x,r)  =^^(X)  < 0 in  we  have  in  this  interval 


i(i(X,r),r,n)  = ^(S^(x),r,n) 


-I 

m=0 


(|S,(X)|r2/2)"'  (/|S,(X)|r)'’ 


(3.3.39) 


m! 


m 


n (2k+l)+2n) 
k=l 


n (2k+l) 
k=0 


and 


Y(S(X.r),r,n)  * Y(S^(x),r.n) 


( i 


m 


ni«0  m 

m!  n (2k-l-2n) 
k=l 


-) 


n (2k-i) 
k=0 


(’^|liU)|r) 


n+i 


(3.3.40) 


The  functions  defined  by  (3.3.39)  and  (3.3.40)  do  not  satisfy  Bessel's 
differential  equation,  but  they  may  be  expressed  in  terms  of  Bessel  functions 
of  a purely  imaginary  argument.  This  is  the  way  they  are  developed  in  our 
computer  program. 

We  begin  our  search  for  eigenvalues  by  finding  the  unique  solution 


Z^(r,X)  of  equation  (3.3.29)  which  satisfies  the  condition 


Lim 


r*0 IT 


= 1 


(3.3.41) 
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and  the  requirement  that  Z^(r,X)  and  K(r)(3/3r)Z^(r,X)  be  contin- 
uous on  [0,  R|^),  We  then  define  a function  of  X by  the  rule, 

Fn(p,c,K,x)  = (3.3.42) 

where  Rj^  is  the  radius  of  the  bounding  sphere  of  the  scattering  body.  We 
then  use  a root-binding  routine  to  find  for  each  n an  ascending  sequence, 

k (3.3.43) 

of  positive  real  numbers  such  that  X = X^^  implies  that 

F„(p,c,K,X)  = 0.  (3.3.44) 

These  numbers  X^^  are  the  eigenvalues  associated  with  th'  >''i4t  transfer 
problem  and  have  the  units  of  reciprocal  seconds.  The  numbers  t^^  = 

1/X(^  k)  estimate  the  time  needed  for  the  (n,k)th  mode  of  the  temperature 
solution  to  decay  to  (1/e)  times  its  original  value. 

3.3.3.  Eigenfunction  Computation.  We  assume  here  that  the  eigenvalue 
X = x^j^  that  we  are  using  to  develop  the  radial  eigenfunction  is  known  and 
use  the  initial  condition  (3.3.26)  and  the  regularity  conditions  (3.3.17)  and 
(3.3.1«)  to  uniquely  determine  the  eigenfunction  Z^^ 

A first  step  in  carrying  this  out  is  the  determination  of  the  eigenfunc- 
tion coefficients  A^.  and  used  in  expressing  the  eigenfunction  Z(r,X) 
by  the  relation  (3.3.25).  We  observe  that  and  are  given  by  equa- 
tions (3.3.26)  and  (3.3.27)  and  that  if  A^  and  B.  are  known,  then 


4 


^4 


:: 

det 

Ai®i[x,Ri,n)  + B^3^(x,R^,n)  »ii) 

Ai^i  (x,Ri  ,n  j + B^B^(x,R^,n)  8^^.^(x,R^  ,n) 

(3.3.45) 

t;  ■■ 

^i+i  ' 

[x ,R^  ,n ) 

i 

det 

B.  = 

Oi^l(x,Ri,n)  A^a^(x,R^,n)  + B^3^(x,R^,n) 

ai+^(x,Ri,n)  A^a.(x,R.,n)  + B^i^(X,R^,n) 

(3.3.46) 

1+1 

(x,R^ ,n) 

where 

1 

o^(X,r,n)  • J(S^(x),r ,n). 

(3.3.47) 

“ •* 

L * . 

o^(x,r,n)  » K^(3/3r)J(S^(x),r,n), 

(3.3.48) 

1 

and 

3^(x,r,n)  « I(l^(x),r,n), 

(3,3.49) 

■ 

B^(x,r,n)  = (3/3r)Y($^ (X),r ,n). 

(3.3.50) 

defines  the  entries  in  the  numerators  of  (3,3.45)  and  (3,3.46)  and  where 


4^(A,R^,n)  = fn)3^^j[A»R^ ,n) 


(3.3.51) 


defines  the  determinant  of  the  matrix  multiplying  the  column  vector  whose 
entries  are  and  8^^^.  Thus,  the  relations  (3.3,45)  and  (3.3.46) 

determine  and  for  all  i e {1,...,N}.  Consequently,  if  X = 
the  eigenfunction  ,^j(r)  has  for  r in  the  explicit  repre- 


sentation 


Z(n,k)(r)  = A^0(S^(X),r,n)  + Y(S^ (X),r,n)  (3.3.52) 
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Y depend  on  whether  or  not 


• .t 


ll  (R._j  < r < R^)  (3.3.53) 

ner  indicated  in  Section  3.2. 


f- 

13 


► *•  * 

• a 


3.4.  Details  of  the  Temperature  Computation 
Including  Complex  Pulse  Heating  Schemes 

jgrlfiS  Expansion  of  the  Temperature.  In  this  section  we  describe 
the  computational  procedure  for  determining  the  solution  u of  equation 
( .3.2)  under  the  assumption  that  we  know  the  eigenvalues  X/  i,\  and  eigen- 
unc  ions  described  in  Section  {3.3.1).  Now  that  this  is  done  we 

express  the  solution  u(r.0,,j,,t)  by  the  series 


u(r,8,<|),t)  = 
n 

I I 

k 


III  3iJ"'’"^(t)P[J(cos(0))cos(m^)Z,  ,^.(r) 
=1  n=0  m=0  ” 


(3.4.1) 


where  a,J"’*''^(t)  is  defined  by  equation  (3.21)  and 


(3.4.2) 


with  the  operators  and  being  defined  by  equations  (3.3.22), 

(3.3.11)  - (3.3.12),  and  (3.3,9)  - (3.3.10)  respectively.  Almost  all  of  the 
computing  time  is  taken  up  in  the  computation  of  the  coefficients  b||’^’'^^(t), 
defined  by  equation  (3.3.23),  that  are  used  in  expanding  the  source  function 
(S/pc).  While  each  of  these  represents  the  result  of  a triple  integration, 
the  total  running  time  is  still  only  between  3 and  4 min  on  an  IBM  360  for 
results  which  are  good  to  within  the  capabilities  of  experimental  measurement. 
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3.4.2.  Complex  Pulse  Heating  Scheme.  We  wish  to  consider  a pulse  heat- 
ing scheme  (e.g..  Figure  3.4.1)  in  which  a group  of  pulses  with  a duty  time 
and  period  followed  by  a period  of  quiescence  defines  a function  that  is 
periodic  with  respect  to  the  total  duty  time  of  the  pulse  group  plus  the 
length  of  time  of  the  quiescent  period. 

More  precisely  the  time  profiles  we  consider  include  time  harmonic  radia- 
tion whose  basic  frequency  is  that  of  a radar  transmitter  multiplied  by  a 
function  of  time  (^cl’^p’^p‘^P’^R*^^  defined  for  0 < < Tp  <.  NpTp  _<  Tp 

and  Tp  > 0 by  the  initialization  rule 


0 
i 

0 
0 

and  the  periodicity  rules, 

if  t + r < Tp  and  t + T„  < N_T  and 

p - R p _ p p 

if  t + Tp  j<  Tp,  where 


t > Tp, 

0 _<  t _<  and  t _<  Tp, 
Th  < t < T„, 


NpTp  < t < Tp, 


(3.4.3) 


(3.4.4) 


(3.4.5) 


T^  = the  pulse  duration, 

Tp  = the  intrapulse  group  period, 

Np  = the  number  of  pulses  per  group. 


51 


OVERfiLL  PICTURE 


AMPLIFIED  PICTURE 


Figure  3.4J.  Conplex  pulse  heating  pattern  typical  of  radar  emissions  with  a 
burst  of  three  pulses  followed  by  a quiet  period  and  with  the 
pattern  being  repeated  periodically.  In  the  above  figure  we 
have  Np  = 3 pulses  per  group,  = .1  millisecond  (ms), 

T = .2  ms,  T = .9  ms,  T„  = 1.6  ms,  and  t = 2.05  ms. 


the  period. 


and 


T 

'R 


the  time  that  the  source  has  been  on. 


t = the  time  of  observation  of  the  radiation  effect. 

The  basic  idea  is  to  assume  that  is  large  enough  that  the  continuous- 
wave  solution  accurately  predicts  the  electromagnetic  field  distribution  and 
consequently  the  source  term  of  the  heat  transfer  equation.  That  is  to  say, 
if 

b['"’"^t)  = (3.4.6) 

then  T = Tj^  implies  that 


3(m,n)(t)  ^ ^m,n) 


r[T/Tp]  N 

l~  I 

k=l  j=l, 


(j-l)Tp  + (k-DTp  + Tj 
exp(-X(t-T))dT  ~ 
(j-l)Tp  + (k-l)Tp 


>^in{[(T  - [T/Tp]Tp)/T  ],N  ) 

- 1“  “ 
j=l 


fftT/Tp]Tp  + (j-DTp  + T^l 
exp(-X(t-T)dT 
[T/Tp]Tp  + (j-l)Tp 

^in(T,{{[T/Tp]Tp  + (min{Np,[(T  - [T/Tp]Tp)/Tp] ) )Tp}  + T^}) 
exp{-X(t-T)dT 

[T/TplTp  (min{Np,t{T  - [T/Tp]Tp)/Tp] ) )Tp 


(3.4.7) 


where  min  is  an  abbreviation  for  the  minimum  function  and  [ ] denotes  the 
greatest  integer  function([xl  denotes  the  largest  integer  not  exceeding  x). 


Some  changes  of  variables  and  introduction  of  notation  will  make  it  easier  to 
develop  computer  code  to  evaluate  the  preceding  three  integrals.  We  therefore 

define 

r(k,j)  = (j-l)Tp  + (k-l)Tp  , 

(3.4.8) 

s(a)  = [T/TplTp  + (M)Tp. 

(3.4.9) 

M(T.Np)  = mn({(T-(T/Tp)Tp)/Tpl,Np), 

(3.4.10) 

l(T.Np.Tp.Tp)  = K(T.Np)Tp  . iT/TpjTp, 

(3.4.11) 

0 if  [(T-[T/TplTp)/Tp]  > Np. 

X = x(Np,T,Tp.Tp)  =1 

— 1 otherwise 

(3.4.12) 

t(T,Np,Tp.Tp)  = min  (T,t(T.Np,Tp,Tp))  + T^jx(Np,T,Tp.Tp), 

(3.4.13) 

Unpj  % 

5(To.Tp,N  ,Td»T)  = I 1 exp(Xr(k,j))  = 

^ ^ k=l  j=l 

exp(XN  T )-l  exp(x[T/Tp]Tp)-l 

{ }{  - - 1. 
exp(^Tp)-l  exp(xTp)-l 

(3.4.14] 

and 

M(T.Np) 

r(T  Tp,N  ,T^,T)  - y exp(Xs(j))  = 

~ j=l 

exp(xM(T,Np]Tp)-l 

exp(x([T/TplTp)){ s , ■■  )• 

(3.4.15 

Putting  all  this  together  we  find  that 


X 

* nTp.VNp.T^.T)) 


exp(xt(T,N  ,T  ,Tp 
+ exp  (-At) ::: 


1 

))  - exp[At(T,Np,Tp,Tp)) 
_ 


This  completes  the  discussion  of  our  temperature  computation  method 


3.5.  Simulated  Biostructures 


3,5.1.  Description  of  Structures  to  be  Studied.  In  a previous  paper 
[2],  the  authors  made  a study  of  the  thermal  response  of  a ball  of  muscle- 
equivalent  material;  this  study  '"s  extended  in  this  pape*'  to  multi-layer 
simulated  biological  structures.  In  [2]  analytical  results  were  compared  with 
measurements  made  with  Vitek-Model  101  Electrorhormia  monitor;  in  this  paper 
we  compare  the  shape  of  the  thermal  response  curve  with  the  electromagnetic 
power  density  curve  which  serves  as  a source  term  for  the  heat  equation.  •'Je 
study  a one-layer  structure  with  blood  flow  at  a higher  frequency  (4.5  GHz) 
than  was  considered  in  [2j,  three-layer  simulated  fetal  structures  with  and 
without  blood  flow,  and  six-layer  "i;.iulated  cranial  structure  with  blood  flow. 

For  one-layer  structures  Figures  2.2.r-2.2.4  show  tho  agreement  between 
theory  and  experimental  measurement;  from  the  results  described  in  Figure 
3.5.1  we  see  that  there  are  striking  resonance  effects  in  a simply  one-layer 
structure  exposed  to  4. 5- GHz  radiation;  we  demonstrate  by  the  results  shown  in 
Figures  3.5.2  and  3.5.3  that  the  temperature  distribution  curve  has  a shape 
very  similar  to  that  of  a power  density  distribution — particularly  when  the 
exposure  time  is  short. 

Next  we  treat  simulated  fetal  structures.  M.  J.  Edwards  [4]  observation 
that  microwave  heating  of  rat  embryos  can  cause  teratogenic  effects  suggests 
that  a quantitative  analysis  of  a simulated  fetal  structure's  response  to 
microwave  radiation  may  assist  in  the  assessment  of  the  potential  hazard  of  a 
source  of  microwave  radiation.  We  use  a three-layer-  model  whose  layers 
consist  of  fetal  tissue,  amniotic  fluid,  and  maternal  tissue  ii  simulating  the 
response  of  the  fetus  to  microwave  radiation.  Figure  3.5.4  shows  the  power 
density  across  the  diameter  of  the  three-layer  structure;  this  diameter  coin- 
cides with  the  z-axis  of  a coordinate  system  whose  origin  i?  the  center  of 


sphere;  the  wave  is  assumed  to  propagate  in  the  direction  of  the  positive 
z-axis.  The  temperature  distributions  across  the  parts  of  the  x,  y,  and  z- 
axes  of  the  structure  within  its  interior  after  a 1-hr  exposure  are  given  in 
Figures  3. 5. 5-3. 5. 7 where  we  include  blood  flow.  Figures  3.5.8-3.5.15  show 
how  this  temperature  distribution  changes  as  exposure  time  increases  when  we 
don't  assume  removal  of  heat  by  blood  flow. 

Finally  we  consider  a six-layer  simulated  cranial  structure  exposed  to 
microwave  radiation.  Figure  3.5.16  shows  the  power  density  across  a six-layer 
simulated  structure  exposed  to  800-MHz  radiation.  Figures  3.5.17  and  3.5.18 
show  the  thermal  response  of  the  structure  to  800-MHz  microwave  radiation  for 
30-s  and  1-min  exposures. 

3.5.2.  Microwave  Heating  of  a Muscle-Equivalent  Sphere.  In  this  section 
we  study  the  manner  in  which  a microwave-induced  temperature  profile  is 
smoothed  as  exposure  time  increases.  We  conclude  that  short-time  temperature 
measurements  would  serve  as  an  adequate  means  of  validating  computer 
predictions  of  internal  field  distributions  even  when  there  are  resonance 
effects  which  cause  the  power  density  profile  (e.g..  Figure  3.5.1)  to  have 
many  relative  maxiinums  and  minimums;  this  particular  assertion  is  valid  for 
continuous-wave  exposure,  but  is  not  established  for  a general  pulse  exposure 
scenario.  The  thermal  response  for  5-s  and  1-min  exposures  is  shown  in 
Figures  3.5.2  and  3.5.3.  The  electromagnetic  field  strength  for  the  results 
portrayed  in  Figures  3.5. 1-3. 5. 3 was  194.09  V/m  or  10  mW/cm^. 

Table  3.5.1  defines  the  parameters  used  in  making  the  computer  runs. 
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TABLE  3.5.1.  PARAMETERS  FOR  ONE-LAYER  MUSCLE  EQUIVALENT 
SPHERE  EXPOSED  TO  4500-MHZ  RADIATION 


ELECTRICAL  PROPERTIES 

Radius  of  bounding  Conductivity 

Tissue  type  sphere  (cm)  Relative  permittivity  (mhos/meter) 


Muscle  3.3 


48.25 


2.75 


THERMAL  PROPERTIES 
(centimeter-gram-second  units) 
Thermal  Specific 

Tissue  type  conductivity  Density  heat 


Blood  flow 
cooling 


Muscle  .00126  1.050  .883  .00186 

3.5.3.  Microwave  Heating  of  a Simulated  Fetal  Structure.  To  estimate 
the  potential  hazard  of  a source  of  microwave  radiation,  we  have  made  a 
simulated  fetal  structure  comprised  of  three  tissue  regions  delimited  by 
concentric  spheres.  The  parameters  used  in  the  computer  runs  are  given  in 
Table  3.5.2.  The  field  strength  used  in  the  runs  was  194.09  v/m,  which  is 
equivalent  to  10  mW/cm^. 


TABLE  3.5.2.  PARAMETERS  DEFINING  A SIMULATED  FETAL 
STRUCTURE  EXPOSED  TO  1000-MHz  RADIATION 


Tissue  type 


ELECTRICAL  PROPERTIES 
Radius 

of  bounding  Relative  Conductivity 

sphere  (cm)  permittivity  (mhos/meter) 


Fetal 

1.6 

50.5 

1.65 

Amniotic  fluid 

2.8 

72.0 

2.00 

Maternal 

3.3 

50.5 

1.65 

THERMAL  PROPERTIES 
(centimeter-gram-second  units) 

Tissue  type 

Thermal  Specinc 

conductivity  Density  heat 

Blood  flow 
cooling 

Fetal 

.00126 

1.050 

.883 

.00186 

Amniotic  fluid 

.00124 

1.007 

.998 

.00000 

Maternal 

.00126 

1.050 

.883 

.00186 
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We  get  some  qualitative  information  (e.g.,  the  shielding  effect  of  the 
amniotic  fluid)  about  the  vulnerability  of  the  fetus  to  microwave  exposure 
from  the  computer  results  for  the  simple  model  given  in  this  paper.  Since 
Edwards  [4]  suggests  that  thermal  pulses  can  affect  cell  cycles  and  that  thei'e 
are  teratogenic  effects  associated  with  elevated  fetal  temperatures,  it  would 
probably  be  valuable  to  carry  out  this  analysis  for  a whole-body  model  and  use 
(1)  the  smallest  fetal  temperature  known  to  cause  abnormal  fetal  development 
and  (2)  the  computer  model  for  predicting  fetal  temperature  excursions  as  a 
definitive  way  of  stating  that  a particular  source  of  microwave  radiation  is  a 
potential  health  hazard. 

Figure  3.5.4  shows  the  power  density  across  a simulated  fetal  structure 
when  the  exposure  was  carried  out  in  the  manner  described  in  Figure  2.1. 
Figures  3. 5. 5-3. 5. 7 show  the  thermal  response  of  the  simulated  fetal  structure 
after  a 1-hr  exposure,  and  Figures  3.5.8-3.5.15  show  how  the  temperature 
distribution  across  the  structure  changes  with  time  when  there  is  no  removal 
of  heat  by  an  autothermal  regulatory  process.  While  our  autothermal 
regulatory  process  model  is  based  on  actual  physiological  parameters  relating 
to  blood  flow,  it  can  at  best  be  considered  phenomenological  since  we  have  not 
modeled  the  details  of  the  flow  of  blood  through  vessels  in  the  tissue  and 
have  in  essence  only  added  a dissipative  term  to  the  heat  equation.  However, 
Figures  3. 1.5-3. 5, 7 suggest  that  there  is  in  the  blood  flow  case  a net  heating 
of  the  amniotic  fluid  due  to  the  absence  of  autothermal  regulation. 
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3.5.4.  Microwave  Heating  of  a Simulated  Cranial  Structure.  In  this 
section  we  study  the  response  of  a simulated  cranial  structure  to  microwave 
radiation.  The  manner  in  v/hich  the  structure  is  exposed  is  described  in 
Figure  2.1.  The  power  density  in  the  simulated  cranial  structure  is  shown  in 
Figure  3.5.16.  The  observed  thermal  response  after  30-s  and  l-min  exposures 
to  800-MHz  radiation  is  described  in  Figures  3.5.17  and  3.5.13. 

The  parameters  used  in  carrying  out  thise  computations  are  described  in 
Table  3.5.3. 


TABLE  3.5.3.  PARAMETERS  DEFINING  A SIX-LAYER  SIMULATED 

CRANIAL  STRUCTURE  EXPOSED  TO  800-MHz  RADIATION 

ELECTRICAL  PROPERTIES 
Radius 

of  bounding  Relative  Conductivity 

Tissue  type  sphere  (cm)  permittivity  (mhos/meter) 


Brain 

2.68 

33.76 

0.960 

CSF 

2.88 

79.47 

1.740 

Dura 

2.93 

45.64 

1.230 

Fat 

3.13 

5.61 

0.096 

Bone 

3,20 

5.61 

0.096 

Skin 

3.30 

45.64 

1.230 

THERMAL 

PROPERTIES 

(centimeter-gram-second 

units) 

Thermal 

Sr  r./ic 

Blood  flow 

Tissue  type 

conductivity 

Density 

heat 

cooling 

Skin 

.0012300 

1.0000 

.900 

.001002242 

Fat 

.0003822 

.8500 

.600 

.000000000 

Bone 

.0027780 

1.5000 

.380 

.000000000 

Dura 

.001230 

1.0000 

.900 

.000000000 

CSF 

.001240 

1.0069 

.998 

4.498x10-® 

Brain 

.001260 

1.0500 

.883 

.00743742 

We  see  from  the  results  described  in  Figures  3.5.17  and  3.5.18  that  there 
is  fairly  rapid  smoothing  of  the  temperature  distributions.  Indeed,  we  have 
assumed  mathematically  that  both  the  temperature  excursion  u and 
K{x,y,z)grad(u),  where  K = K(x,y,2)  is  the  thermal  conductivity,  are  contin- 
uous. Thus,  since  in  our  calculation  K is  constant,  we  see  that  u and  its 
derivatives  are  continuous.  The  electrical  properties  of  the  six  tissue  types 
are  given  in  [7].  The  model  is  capable  of  predicting  the  microwave-inducei. 
temperature  excursion  when  there  is  blood  flow  in  some  of  the  layers,  and 
typical  values  of  these  blood  flow  parameters  can  be  obtained  from  [6]. 
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igure  3.5 


2.  Thermal  response  of  a muscle-equivalent  sphere  to  a 1-min  expo- 
sure to  4.5-GHz  continuous-wave  radiation  with  a pov/er  of  10 
mW/cm^,  Parameters  defining  the  problem  are  given  in  Table 


teuperrture  rise  coeg 
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Figure  3.5.3.  Thermal  response  of  a muscle-equivalent  sphere  to  a 5-s 
exposure  of  4.5-GHz  continuous-wave  radiation  with  a power  cf 
10  mW/cm^.  Parameters  defining  the  problem  are  oiven  in  Table 
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Figure  3.5,4.  Power  density  across  the  z-axis  of  a simulated  fetal  structure 
exposed  to  1-GHz  continuous-wave  microwave  radiation  with  a 
power  of  10  mW/cm^.  Parameters  defining  the  problem  are  given 


in  Table  3.6.2. 
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Figure  3.5 


5.  Thermal  response  of  a simulated  fetal  structure  to  a 1-hr 
exposure  to  1-GHz  radiation  with  a power  of  10  mW/cm^.  The 
temperature  is  computed  across  the  x-axis.  The  orientation  of 
the  axes  is  given  in  Figure  2.1.  ‘’'he  parameters  defining  the 

problem  are  given  in  Table  3.5.2. 
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Figure  3.5 


9,  Temperature  distribution  along  the  z-axis  of  a simulated  fetal 
structure  exposed  to  1-GHz  (10  mW/cm^)  radiation  for  1 min. 
The  parameters  defining  the  problem  are  given  in  Table  3.5.2, 
except  that  all  blood  flow  terms  are  set  to  zero. 
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Figure  3.5 


10.  Temperature  rise  along  the  z-axis  of  a simulated  fetal  struc- 
ture exposed  to  1-GHz  (10  mW/cm^)  radiation  for  15  min.  The 
parameters  defining  the  problem  are  given  in  Table  3.5.2, 
except  that  all  blood  flow  terms  are  set  to  zero. 
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Figure  3.5.11.  Temperature  rise  along  the  z-axis  of  a simulated  fetal  struc- 
ture exposed  to  1-GHz  (10  mW/cm^)  radiation  for  1 hr.  The 
parameters  defining  this  problem  are  given  in  Table  3.5.2, 
except  that  all  blood  flow  terms  are  set  to  zero. 
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Figure  3 


5.12,  Temperature  rise  along  the  z-axis  of  a simulated  fetal  struc- 
ture exposed  to  1-GHz  (10  mW/cm^)  radiation  for  2 hr.  Th( 
parameters  defining  this  problem  are  given  in  Table  3,5.2 
except  that  all  blood  flow  terms  are  set  to  zero. 
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Figure  3.5.13.  Temperature  rise  along  the  .T-axis  of  a simulated  fetal  struc 
ture  exposed  to  l-6Hz  (10  mW/cm^)  radiation  for  3 hr.  Th 
parameters  defining  this  problem  are  given  in  Table  3.5.2 
except  that  all  blood  flow  terms  are  set  to  zero. 
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Figure  3.5.14.  Temperature  rise  along  the  z-axis  of  a simulated  fetal  struc- 
ture exposed  to  1-GHz  (10  mW/cm^)  radiation  for  4 hr.  The 
parameters  defining  this  problem  are  given  in  Table  3.5.2, 
except  that  all  blood  flow  terms  are  set  to  zero. 
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Figure  3.5 


16.  Power  density  along  the  z-axis  of  a six-layer  simulated  cranial 
structure  exposed  to  800-MHz  radiation  with  a power  of  10 
mW/cm^.  The  parameters  defining  this  problem  are  given  in 


Table  3.5.3. 
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Figure  3.5.17.  Thermal  response  of  a six-layer  simulated  cranial  structure 
exposed  tc  800-MHz  radiation  for  3 min.  The  parainecern  defin- 
i.g  this  problem  are  given  in  Table  3.5.3. 
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4.  PROGRAM  DESCRIPTION 


4.1.  Purpose  of  the  Program 

The  computer  program  described  in  this  report  will  predict  the  thermal 
response  of  an  autothermal ly  regulated,  spherically  symmetric,  dielectric  body 
with  a finite  microwave  conductivity  to  a time-harmonic  source  of  microwave 
radiation.  The  calculation  can  be  carried  out  at  any  point  in  the  interior  of 
this  body  at  any  positive  time.  We  also  allow  the  source  to  be  pulsed  in  the 
sense  that  the  source  of  time  harmonic  radiation  may  be  turned  on  and  off  in  a 
complex  manner  such  as  that  described  in  Figure  3.4.1.  The  scattering  body 
consists  of  from  one  to  six  homogeneous  regions  bounded  on  the  outside  by  a 
sphere  of  finite  but  positive  radius;  a description  of  a six-layer  body  is 
shown  in  Figure  3.5.1.  The  radiation  source  is  given  by  an  amplitude  Eg  in 
volts  per  meter  and  a frequency  FREQ  in  megahertz.  The  calculation  is  carried 
out  by  the  evaluation  of  an  analytical  expression  involving  an  infinite  linear 
combination  of  spherical  harmonics.  The  coefficients  in  this  infinite  linear 
combination  are  functions  of  time  computed  from  a set  of  eigenvalues  and  a 
kn'iwlecge  of  the  manner  in  which  the  microwave  source  has  been  turned  on  and 
off.  We  also  permit  a nonzero  heat  removal  term  BFRP  in  one  or  more  of  the 
layers. 


4.2.  Accessing  the  Program  from  the  Library 
The  Job  Control  Language  needed  to  access  the  program  from  the  library  is 
shown  in  Figure  4.2.1.  The  data  deck,  whose  preparation  will  be  explained  in 
Section  4.2,  goes  between  the  card  marked 
//GO.SYSIN  DD  * 
and  the  card 
/* 

The  space  requirement  specified  by  the  expression 
REGION. G0=252K 

will  not  change  out  if  one  desires  temperature  information  at  more  points  one 
needs  to  increase  the  running  tiine  parameter, 

TIME.G0=4 

and  the  effective  time  parameter, 

EFFTIME=10 

which  represents  elapsed  time  in  our  timesharing  computing  system.  The  para- 
meter values  listed  in  Figure  4.2.1  were  adequate  to  determine  the  temperature 
excursions  for  a single  fixed  time  at  13  different  locations  in  space. 

//HBR16TKS  J0B(3H01,B02O,,,,,,,00), 'HBM3790(310R  COHOON ' ,CLASS=C , 

//  MSGLEVEL=(2,0) 

/♦JORPARM  RESTART, SINGLE,  EFFTIME=10,L INES=1 ,CARnS-n 
//STEPl  EXEC  PL0TG0,PR0GRAM=HBR16TRF, REGION. G0=252K, TIME. Gn=4 
//GO.SYSIN  DD  * 

DATA  CAPOS  GO  HERE  **♦ 

/* 

Figure  4.2.1.  Job  control  language  for  calling  the  microwave  thermal  response 

program  from  the  library. 
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4.3.  Glossary  of  Variables  and  Their  Meaning 


All  FORTRAN  variables  used  for  input  and  output  and  important  internal 
FORTRAN  variables  are  listed  here  in  alphabetical  order.  We,  with  each  vari- 
able, give  an  explanation  of  its  meaning.  These  variables  are: 

ALPNP(NNN)  = the  ALPHA  SUB  (L,P)  coefficient  used  in  formula 
(3.2.2)  to  expand  the  electromagnetic  field  with 
NNN  = (NREG-1)*NMIN+NN 


and 


ALPNP(NNN)  = ALPHA  SUB  {NN,NREG)  where  NN  is  the  spherical 
Bessel  function  order  and  NREG  is  the  layer  number. 

ANP  = the  A SUB  (L,P)  coefficient  used  in  expanding  the  electric  field  and 
which  appears  in  formula  (3.2.2)  which  is  stored  in  a single  array  with 
ANP((NREG-1)*NMIN+NN)  corresponding  to  the  coefficient  A SUB  (NN,NREG)  where 
NN  indexes  the  spherical  Bessel  functions  used  in  the  expansion. 

BETNP(NNN)  = the  BETA  SUB  (L,P)  coefficient  used  in  formula  (3.2.2)  to  expand 
the  electromagnetic  field  with 
f NNN  = (NREG-1)*NMIN+NN 

and 

RETNP(r;NN)  = BETA  SUB  (NN,NRE6) 

where  NN  is  the  spherical  Bessel  function  order  and  NREG  is  the  layer  number. 

BFRP(I)  -•  the  blood  flow  radial  perviousness  term  ot'  the  number  of  grams  of 
blood  per  gram  of  tissue  per  second,  with  a typical  value  for  brain  tissue 
being  .0122 

BNP  = the  B SUB  (L,P)  coefficient  used  in  formula  (3.2.2)  to  expand  the 
electromagnetic  field  and  which  is  stored  in  a single  array  with 
BNP( (NRTG-1 )*NMIN+NN)  corresponding  to  the  coefficient  B SUB  (NN,NREG)  where 
NN  indexes  the  spherical  Bessel  functions  used  in  the  expansion. 

RP(1)  = the  product  of  BFRP(I),  the  number  of  grams  of  blood  per  gram  of 
tissue  per  second,  CBP  = .98  = the  specific  heat  of  blood  in  calories  per  gram 
degree  Centigrade,  and  the  density  RHOBP  = 1.06  = the  number  of  grams  of 
tissue  per  cu'jic  centimeter  of  tissue. 

CALL  BJYH(BJNP,BHNP,0,NC ,ST0PR,MAX)  = a call  to  a subroutine  which  determines 
the  values  of  spherical  Bessel  fi-nctions  of  the  first  kind  BJNP  and  spherical 
Hankel  functions  BHNP  at  the  complex  argument  0.  We  attempt  to  generate  up  to 
MAX  such  functions  as  we  are  limited  by  STOPR  and  we  end  up  putting  only  NC 
such  functions  in  the  array. 

CALL  COEF  = a call  to  a subroutine  which  produces  the  expansion  coefficients 
A SUB(L,P),  alpha  SUB(L,P), 

B SUB(L,P)  and  BETA  SUB(L,P) 

used  in  expanding  the  electromagnetic  field  using  equation  (3.2.2) 


mratfrUrfi-.  n 
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CALL  DRTMI{X,F,FNCAL,SL,SR,W,V,E,NITR)  = the  call  to  the  bisection  routine 
ORTMI  which  returns  a value  X such  that  FNCAL(X)  = F =0  to  within  an  accuracy 
of  E with  less  than  NITR  iterations  where  FNCAL(SL)=W,  FNCAL(SR)=V  and  W and  V 
are  on  opposite  sides  of  0 on  the  real  line. 

CALL  EPROP(FREQ,ITIS(I),EP,SIG)  = a call  to  a subroutine  which  determines  the 
relative  permittivity  EP  and  microwave  conductivity  SIG  of  tissue  type  ITIS(I) 
at  the  microwave  frequency  FREQ. 

CALL  PL  = a call  to  a subroutine  which  computes  the  array  P of  associated 
Legendre  polynomials  of  the  first  kind  and  order  1 and  an  array  DP  of  their 
derivatives. 

CALL  TERM(NCK,T,KEY)  = a call  to  a subroutine  which  computes  the  I**L  multi- 
plied by  T appearing  in  formula  (3.2.2)  based  on  its  preceding  value,  where 
the  value  of  NCK  ranges  from  1 to  4 since  1**1,  1**2,  1**3,  1**4  ranges  over 
all  possible  values  of  the  square  root  of  (-1)  raised  to  a power,  and  where 
KEY  takes  on  the  value  1 or  0 depending  on  where  in  the  process  of  summing  the 
series  we  are  computing  I**L  multiplied  by  the  complex  term  T appearing  in 
equation  (3.2.2). 

CP(I)  = the  tissue  specific  heat  in  c'lories  per  gram  per  degree  centigrade. 

DEN(f>'SBF,M2,NRT)  = the  integral  of  the  square  of  the  radial  eigenfunction 
multiplied  by  the  square  of  the  radial  coordinate,  the  density,  and  the  spe- 
cific heat  from  zero  to  the  outer  radius  of  the  scatterer. 

EO  = the  strength  of  the  incident  electric  field  vector  which  may  be  read  in 
a certain  number  of  milliwatts  per  square  centimeter  if  lEO  = 1 but  must  be 
expressed  in  volts  per  meter  if  lEO  « 0,  where  we  understand  that  if  lEO  = 1, 
then  EO  will  be  converted  internally  into  volts  per  meter. 

EPHI  = the  phi  component  of  the  electric  field  vector  when  the  electric  field 
is  expressed  in  spherical  coordinates  and  which  corssequently  represents  a 
tangential  field  component  when  the  point  at  which  the  field  is  being  computed 
is  on  a sphere  defining  a boundary  of  the  body  being  heated  by  microwaves. 

EPS  = the  relative  error  associated  with  the  expansion  of  the  electromagnetic 
field. 

EPSP(I)  = the  relative  dielectric  constant  of  the  Ith  tissue  layer  at  the 
frequency  FREQ  of  the  incoming  radiation. 

ERAO  = the  radial  component  of  the  electric  vector  in  volts  per  meter  where  we 
assume  that  we  have  expressed  the  field  vectors  in  the  spherical  coordinate 
system  and  which  consequently  represents  the  component  of  electric  vector  that 
is  perpendicular  to  a boundary  layer  when  the  point  at  which  the  electric 
field  is  being  computed  is  on  a sphere  defining  a boundary  of  the  body  being 
heated. 

ETHETA  = the  theta  r''-.ponent  of  the  electric  field  vector  v/hen  the  electric 
field  is  expressed  in  spherical  coordinates  and  which  consequently  represents 
a tangential  field  component  when  the  point  at  which  the  field  is  being 
computed  is  on  a sphere  defining  a boundary  of  the  body  being  heated  by  micro- 
waves. 
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ETIME(NRT)  = the  time  profile  function,  defined  by  dividing  the  right  side  of 

equation  (3.4.7)  by  b-SUB-BAR-SUB-K-SUB-(M,N)  or  which  describes  the 
radar  pulse  emission  patterns. 

F = the  factor  in  front  of  the  integral  on  the  right  side  of  equation  (3.3.11) 
in  general  equal  to  (2n+l ) ( (n-m) ! )/( (n+m) ! ) multiplied  by  the  factor  in  front 
of  the  integrals  on  the  right  side  of  equation  (3.3.9)  or  (3.3.10)  whichever 
is  appropriate. 

FKP(NREG)  = the  complex  electromagnetic  propagation  constant  associated  with 
layer  NREG  which  is  defined  by  equation  (3.2.7). 

FREQ  = the  frequency  of  the  incoming  radiation  in  megahertz  or  millions  of 
cycles  per  second. 

FUNCTION  ALP(N,M,X)  = a function  subroutine  computing  the  associated  Legendre 
function  of  the  first  kind  of  degree  N and  order  M at  the  point  X with  the 
restriction  that  N and  M are  nonnegative  integers  and  M does  not  exceed  N. 

FUNCTION  FNCAL(EIGV)  = a function  subroutine  whose  output  is  the  value  of  the 
Newton  cooling  function  defined  by  equation  (3.3.42)  when  LAMDA  = EIGV. 

lEO  = a parameter  for  determining  the  way  the  input  data  EO  is  interpreted 
with  lEO  = 0 meaning  that  EO  is  a certain  number  of  volts  per  meter  and  lEO  = 

I meaning  that  EO  is  a certain  number  of  milliwatts  per  square  centimeter. 

II  = in  the  last  print  statement  an  index  describing  the  number  of  the  data 
card  containing  the  point  at  which  the  temperature  is  being  computed  with  II  = 
1 for  the  point  on  the  first  card  and  II  = NOCR  for  the  point  on  the  last 
car^.. 

ISAR  = a parameter  determining  the  way  that  the  output  data  is  expressed  with 
ISAR  = 0 if  the  predicted  power  density  that  is  printed  next  to  the  predicted 
temperature  is  to  be  expressed  in  milliwatts  per  kilogram,  and  ISAR  = 1 if  it 
is  to  be  expressed  in  watts  per  cubic  meter. 

ITIS(I)  = the  tissue  type  of  the  Ith  tissue  layer  equal  to  1,2, 3, 4, 5, 6,  or  7 
if  the  tissue  type  is  (i)  cerebrospinal  fluid,  (ii)  blood,  (iii)  muscle,  (iv) 
skin  or  dura,  (v)  brain,  (vi ) fat  or  bone,  or  (vii)  yellow  bone  marrow, 
respectively. 

KMAX  = the  number  of  radial  eigenfunctions  associated  with  a given  order  of 
Bessel  function  with  the  greatest  accuracy  being  achieved  by  setting  KMAX 
equal  to  its  maximum  value  of  25. 

MP  = the  number  of  points  to  be  used  in  the  Gauss  quadrature  integration 
scheme  that  executes  the  radial  transform  defined  by  equation  (3.3.22)  with 
this  number  being  one  of  32,  48,  64,  or  80  and  with  the  larger  numbers  giving 
the  more  accurate  results. 

MPl  = the  number  of  points  used  in  the  Gauss  quadrature  scheme  that  performs 
the  Legendre  transform  defined  by  equation  (3.3.1)  with  the  number  being  32 
Of  48  and  where  the  latter  number  gives  the  most  accurate  results. 
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NC  = the  maximum  number  of  Bessel  functions  available  based  on  the  value  of 
the  particular  point  at  which  the  field  is  being  computed  and  the  microwave 
electrical  properties  of  the  layer  in  which  the  point  is  located. 

NMAX  = the  number  of  orders  of  spherical  Bessel  functions  that  will  be  used  to 
describe  the  raoial  variation  of  the  microwave  radiation-induced  temperature 
excursion  with  the  greatest  accuracy  being  achieved  by  setting  NMAX  equal  to 
its  maximum  value  of  12. 

NMIN  = the  number  of  expansion  coefficients  available  based  on  the  radii  of 
the  spheres  bounding  the  tissue  layers  and  the  microwave  electrical  properties 
of  the  material  in  these  layers. 

NNN  = (NREG  - 1)*NMIN  + NN,  where  NN  denotes  the  spherical  Bessel  function 
order. 

NREG  = the  number  of  the  layer  in  which  the  point  at  which  the  temperature  is 
being  computed  is  found. 

NOCR  = the  number  of  spatial  points  at  which  the  input  data  set  is  to  be 
computed,  the  maximum  value  of  the  index  II  of  the  output  temperature  data  for 
a particular  exposure  time,  and  the  number  of  cards  in  the  fourth  input  data 
set. 

NORG  = the  number  of  layers  in  the  model  where  NORG  is  1 if  the  scatterer  is 
a homogeneous  ball  and  where  NORG  equals  its  maximum  value  of  6 if  the  body  in 
which  the  microwave-induced  temperature  is  being  predicted  is  a ball  sur- 
rounded by  five  outer  layers. 

NPniNT(I)  = the  Ith  entry  of  a 5-element  array  containing  allowable  numbers  of 
points  that  may  be  used  in  a Gauss  quadrature  scheme  for  evaluating  expansion 
coefficients. 

NPUL  = the  number  of  pulses  in  a group,  where  for  example  NPUL  = 3 if  the 
radar  emission  pattern  being  modeled  consists  of  3 bursts  of  radiation  fol- 
lowed by  a quiet  period,  3 bursts  and  a quiet  period,  et  cetera. 

NSBF  = FORTRAN  index  equal  to  one  plus  the  order  of  the  Bessel  function  being 
considered  in  the  computation  of  the  microwave-induced  temperature. 

PAV6  = the  total  absorbed  power  di video  by  the  total  volume  of  the  region  in 
which  the  temperature  increase  is  being  predicted  expressed  in  watts  per  cubic 
meter. 

PAVGl  = the  total  absorbed  power  divided  by  the  total  mass  in  kilograms  of  the 
body  in  which  the  temperature  excursion  is  being  predicted  expressed  in  milli- 
watts per  kilogram. 

PCEBF  = the  relative  error  in  temperature  computation  associated  with  using 
one  less  order  of  Bessel  function  but  keeping  the  same  number  of  eigenvalues 
for  each  Bessel  function  order  which,  for  example,  would  mean  computing  the 
temperature  SBFMl  using  NMAX-1  Bessel  functions  and  the  full  KMAX  eigenvalues 
per  Bessel  function  and  defining  PCEBF  = (TRM-SBFMl )/TRM  to  be  this  relative 
error. 
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PCER  = the  relative  error  associated  with  leaving  off  the  last  eigenfunction 
or,  for  example,  using  24  eigenfunctions  instead  of  25  eigenfunctions  for  each 
Bessel  function  order. 

PD  = the  value  of  the  divergence  of  the  Poynting  vector  at  the  point  whose 
spherical  coordinates  are  (SAVER, THETAD,PHID)  which  value  represents  the  num- 
ber of  milliwatts  of  power  being  deposited  per  cubic  centimeter  of  tissue. 

PHID  = the  phi  coordinate  of  the  point  at  which  the  microwave-induced  tempera- 
ture is  to  be  computed,  where  phi  is  the  spherical  coordinate  that  ranges 
between  zero  and  360  degrees. 

R = the  radial  coordinate  of  the  point  at  which  the  temperature  is  being 
computed. 

RHOP(I)  = the  tissue  density  of  the  Ith  tissue  layer  in  grams  per  cubic  meter 
where  typically  RH0P(I)=1. 

SBDP(I)  = the  radius  in  centimeters  of  the  smallest  sphere  containing  the  Ith 
tissue  layer  where  I ranges  from  1 to  N0R6. 

SBEMl  = the  predicted  temperature  obtained  by  using  KMAX  roots  per  Bessel 
function  order  but  only  NMAX-1  Bessel  function  orders  in  approximating  the 
infinite  sum  of  equation  (3.4.1). 

SIGP(I)  = the  conductivity  in  mhos  per  meter  of  the  Ith  tissue  layer  where  I 
ranges  from  1 to  NORG. 

SRMl  = the  estimated  temperatu/e  using  NMAX  Bessel  function  orders  and  KMAX-1 
eigenvalues  per  Bessel  function  order. 

STOPR  = a termination  indicator  for  stopping  the  generation  of  Bessel  func- 
tions based  on  the  fact  that  STOPR  exceeds  the  maximum  absolute  value  of  any 
of  the  spherical  Bessel  functions  of  the  second  kind  used  in  describing  the 
dependence  of  the  induced  and  scattered  electromagnetic  fields  on  the  radial 
variable  with  a typical  value  being  1.E35. 

TBPER  = the  period  of  the  pulse  group  envelope,  where,  for  example,  if  there 
is  a radar  emission  pattern  consisting  of  a burst  of  three  pulses  of  duration 
3*TPER  followed  by  a quiet  period,  a burst  of  three  pulses  followed  by  a quiet 
period,  et  cetera,  then  TBPER  is  equal  to  3*TPER  plus  the  major  quiet  period, 
where,  of  course,  we  defin-^  the  major  quiet  period  to  be  total  quiet  period 
minus  the  time  between  the  individual  pulses  in  the  group  or  as  the  T-SUB-p  in 
Figure  4.3.1. 

TCP(I)  = the  thermal  conductivity  in  calories  per  centimeter  per  degree 
centigrade  per  second  of  the  Ith  tissue  layer  where  I ranges  from  1 to  NCRG 
with  a typical  value  being  .0012. 

TOUT  = the  time  at  i.hich  the  source  of  pulsed  microwave  radiation  is  shut  down 
or  the  T-SUB-R  of  Figure  4.4.1. 
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TDUR  = the  up  time  in  seconds  of  an  isolated  pulse  in  a pulse  group  or  the 
value  of  T-SUB-d  in  Figure  4.4.1. 

THETAO  = the  theta  coordinate  of  the  point  at  which  the  temperature  is  being 
computed,  where  this  is  the  spherical  coordinate  that  ranges  between  zero  and 
180  degrees. 

TIME  = the  time  in  seconds  at  which  the  microwave-induced  temperature  is  to  be 
computed. 

TOTPOW  = the  total  absorbed  power  in  watts  determined  by  carrying  out  an 
energy  balance  on  the  surface  of  the  scatterer  using  the  Poynting  vector  for 
the  incident  and  reflected  radiation. 

TRM  = the  microwave-induceo  temperature  obtained  by  adding  up  terms  in  an 
eigenfunction  expansion  at  the  point  whose  spherical  coordinates  are  specified 
by  the  three-tuple,  (^^AVER,  THETAD,PHID). 

U(NSBF,M2,K)  = the  exparsion  coefficient  which  is  defined  by  equation  (3.4.16) 
at  the  observation  time  TIME  and  which  is  used  in  equation  (3.4.1),  where 
NSBF  is  one  plus  the  order  of  the  Bessel  function,  M2  is  1 or  2 depending  on 
whether  the  ind-ax  of  the  cosine  transform  defined  by  equations  (3.3.9)  and 
(3.3.10)  is  0 or  2,  and  K is  the  index  of  the  eigenvalue  associated  with  a 
given  Bessel  function  order. 

VOL  = the  volume  of  the  body  in  which  the  microwave-induced  temperature  is 
being  predicted  in  cubic  meters. 

XLAMOA(K,NSBF)  « an  element  of  a KMAX  by  NMAX  array  (dimensioned  as  25  by  12) 
which  represents  the  Kth  eigenvalue  associated  with  the  Bessel  function  of 
order  NSFF,  where  each  of  these  numbers  is  used  to  define  a combination  of 
spherical  Bessel  functions  satisfying  the  Newton  cooling  law  at  the  outer 

boundary  with  this  combination  of  Bessel  functions  being  the  eigenfunction 
used  in  the  eigenfunction  expansion  of  the  microwave-induced  temperature. 

XMASS  = the  mass  of  the  scattering  body  in  kilograms. 

XNUM(NSBF,M2,NRT)  = the  transform  of  the  source  term  with  respect  to  its 

spatial  variables. 

ZLAB(l)  = the  first  element  of  an  alpha  array  containing  the  expression 

'W/M**3'. 

ZLAB(2)  = the  second  element  of  an  alpha  array  containing  the  expression 

'MW/KG'. 


88 


4.4.  Input  Data  Preparation 


The  purpose  of  this  section  is  to  tell  a user  how  to  prepare  data  to  run 
the  computer  program  to  predict  the  thermal  response  of  a spherically  symmet- 
ric penetrable  body  to  microwave  radiation.  The  incoming  radiation,  the 
precision  with  which  the  response  to  this  radiation  will  be  calculated,  the 
temporal  envelope  of  the  incoming  radiation  and  the  time  at  which  the  tempera- 
ture response  is  to  be  computed,  and  the  thermal  and  electrical  properties  of 
the  body  in  which  the  temperature  excursion  is  being  predicted  are  described 
in  the  first  three  data  sets.  Data  set  three  is  a multi-card  set  with  the 
number  of  cards  being  equal  to  the  number  of  layers  the  scattering  body. 
The  fourth  data  set  is  the  collection  of  points  at  which  one  seeks  to  compute 
the  temperature;  each  point  is  on  a separate  card. 

The  following  paragraphs  give  details  concerning  the  composition  of  the 
four  data  sets  used  by  the  computer  program.  Figures  at  the  end  give  some 
data  sets  that  direct  th^  program  to  predict  the  microwave-radiation-induced 
temperature  on  the  x,  y,  and  z axes  of  the  sphere. 

Data  set  one  consists  of  a single  card  containing  FREQ,  EO,  STOPR,  NORG, 
NMAX,  KMAX,  MP,  MPl,  lEO,  and  ISAR  which  is  read  in  via  the  statements: 

READ  5, FREQ, EO, STOPR, NORG„NMAX,KMAX,MP, MPl, lEO, ISAR 
5 F0RMAT(30i0. 0,715) 

In  the  above 

FREQ  = the  frequency  of  the  incoming  radiation  in  megahertz, 

EO  = the  strength  of  the  incoming  E-field  in  volts  per  meter  (if  lEO  = 0)  and 
in  milliwatts  per  square  centimeter  (if  lEO  = 1), 

STOPR  = a termination  indicator  for  stopping  the  generation  of  Bessel 
functions  based  on  the  fact  that  STOPR  exceeds  the  absolute  value  of  any  of 
the  spherical  Bessel  functions  Y that  will  be  used  in  describing  the 
dependence  of  the  induced  and  scattered  fields  on  the  radial  variable  with  a 
typical  value  being  1.E35. 


NORG  = the  number  of  layers  in  the  model  where  NORG  is  1 if  the  scatterer  is  a 
homogeneous  ball  and  NORG  equals  its  maximum  value  of  6 for  a ball  surrounded 
by  5 outer  layers, 

NMAX  = the  number  of  orders  of  spherical  Bessel  functions  that  will  be  used  to 
help  describe  the  microwave -radiation-induced  temperature  with  the  greatest 
accuracy  being  achieved  by  setting  NMAX  equal  to  its  maximum  value  of  12, 

KMAX  = the  number  of  radial  eigenfunctions  associated  with  a given  order  of 
Bessel  function  where  the  greatest  accuracy  is  achieved  by  setting  KMAX  equal 
to  its  maximum  value  of  25, 

MP  = the  number  of  points  used  in  the  Gauss  quadrature  scheme  that  carries  out 
the  radial  transform  defined  by  the  formula{3.22)  with  this  number  being  one 
of  32,  48,  64,  or  80  and  with  the  larger  numbers  giving  the  more  accurate 
results, 

MPl  = the  number  of  points  used  in  the  Gauss  quadrature  scheme  that  carries 
out  the  Legendre  transform  defined  by  equation(3.11 ) with  this  number  being  32 
or  48, 

lEO  = a parameter  for  determining  the  way  the  input  data  EO  is  interpreted 
with  lEO  = 0 meaning  that  EO  is  a certain  number  of  volts  per  meter  and  lEO  = 
1 meaning  that  EO  is  a certain  number  of  milliwatts  per  square  centimeter, 

and 

ISAR  = a parameter  determining  the  way  that  the  output  data  is  expressed  with 
ISAR  = 0 if  the  predicted  power  density  is  to  be  expressed  in  milliwatts  per 
kilogram  and  ISAR  * 1 if  the  predicted  power  density  is  to  be  written  out  and 
labeled  as  a certain  number  of  watts  per  cubic  meter. 

Data  set  two  consists  also  of  a single  card  containing  TDUR,  TPER,  TBPER,  TOUT, 
TIME,  NPUl.,  and  NOCR,  which  is  read  in  through  the  statements: 

10  REA0(5, 15, END=350)TDUR, TPER ,TBPER,TCUT, TIME, NPUL, NOCR, IPL1,IPL2 
15  F0RMAT(5D10. 0,215) 

In  the  above 

TDUR  = the  duration  of  the  pulse  or  the  value  of  T-sub-d  in  Figure  4.4.1, 

TPER  = the  .''>riod  in  the  primary  pulse  group  or  the  value  of  T-sut-p  in 

Figure  4.4.1, 

TBPER  = the  period  of  the  pulse  group  envelope  or  the  T-sub-P  in  Figure  4.4.1, 
possible  time  envelope  function  for  incoming  radiation. 

This  is  similar  to  some  radar  emission  patterns. 

TCUT  = the  time  at  which  the  source  is  shut  down  or  the  T-sub-R  in  Figure 

4.4.1, 

TIME  = the  time  at  which  the  microwave-induced  temperature  is  to  be  computed. 
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NPUL  = the  number  of  pulses  per-  group,  where  we  note  that  NPUL=2  in  Figure 
A. 4.1  and  NPUL=3  in  Figure  3.4.1, 

NOCR  = the  number  of  spatial  points  at  which  the  temperature  is  to  be  com- 
puted. 

IPLl  - an  integer  ranging  from  0 to  7,  which  will  indicate  which  of  certain 
plots  of  temperature  across  the  sphere  diameters  coinciding  with  the 
coordinate  axes  will  be  given.  A value  of  IPLl  equal  to 

0 means  no  axis  plots  of  temperature  will  be  produced, 

1 means  a plot  of  temperature  across  the  z-axis  will  be  given, 

2 means  a plot  of  temperature  across  the  x-axis  will  be  given, 

3 means  a plot  of  temperature  across  the  y-axis  will  be  given, 

4 means  combined  results  of  1 and  2 are  given, 

5 means  combined  results  of  1 and  2 are  produced, 

6 means  combined  results  of  2 ano  3 are  produced,  and 

7 means  combined  results  of  1,  2,  and  3 are  given. 


IPL2  = an  integer  ranging  from  0 to  7,  which  will  indicate  which  of  v.ertain 
contour  plots  of  isotherms  will  be  produced  on  the  plotting  file  F0R008.DAT. 
In  the  following  description  the  x-z  plane  refers  to  the  intersection  of  the 
plane  containing  the  x-axis  and  the  2-axis  with  the  interior  of  the  bounding 
sphere.  The  y-z  plane  will  mean  the  plane  containing  the  y and  z axes  or  the 
plane  x * 0.  The  x-y  plane  means  the  plane  z = 0.  A value  of  IPL2  equal  to 

0 means  no  axis  plots  of  temperature  will  be  produced, 

1 means  a contour  plot  in  the  x-z  plane  will  be  given, 

2 means  a contour  plot  in  the  y-z  plane  will  be  given, 

3 means  a contour  plot  in  the  x-y  plane  will  be  given, 

4 gives  the  combined  results  of  1 and  2, 

5 gives  the  combined  results  of  1 and  3, 

6 gives  the  combined  results  of  2 and  3,  and 

7 gives  the  combined  results  of  1,  2,  and  3. 
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igure  4,4,1.  Typical  time  envelope  function  dascribing  some  radar  emission 
patterns.  In  the  above  figure  we  have  Np  = 2 pulses  per 

group-,  Tj  = .025  milliseconds  (ms),  Tp  = .050  ms. 

To  = 


9 ms.  and  T-  = .4125  ms. 


Data  set  three  consists  of  NORG  cards  indexed  by  the  parameter  I 
augmented  from  1 to  NORG  in  a 00  LOOP.  The  Ith  card  contains  SBDP(I), 
EPSP(I),  SIGP(I).  TCP(I),  RHOP(I),  CP(I),  BFRP(I).  and  ITIS(I)  read  in  by 
means  of  the  statements: 


30  F0RMM{7K10.0,I5) 

DO  6£  I = 1,N0RG 

READ  30.  SBOP(I),  EPSP(I),  SFGP{I),  TCP(I),  RHOP(I), 
ICP(I),  BFRP(I),  ITIS(I) 


s & • • 


65  PRINT  70.  I,  SBDP{I),  EPSP(I),  SIGP(I). 

ITCP(I),  RHOP(I),  cP(I).  BFRP(I). 

2TISSUE(ITIS(I)). 

70  F0RMAT(I4,F12.2,F12.2,F13.3.F15.6,F13.4, 

1F10.3,F12.5,3X,A8) 

In  the  above,  the  properties  of  the  Ith  layer  are  specified  by  letting 
SBDP(I)  = its  outer  radius  in  centimeters, 

EPSP(I)  = its  relative  permittivity, 

SIGP(I)  = its  conductivity  in  mhos  per  meter, 

TCP(I)  = its  thermal  conductivity  in  calories  per  centimeter  per  degree  centi- 
grade per  second (typically  TCP (I)  = .0012), 

RHOP(I)  = tissue  density  in  grams  per  cubic  centimeter(typically  RHOP(I)  = 

1). 

CP(I)  = tissue  specific  heat  in  calories  per  gram  degree  centigrade{typically 
CP{I)  = .84), 

BFRP(I)  = the  blood  flow  term  that  is  equal  to  the  product  of  the  number  of 
grams  of  blood  per  gram  of  tissue  per  second,  the  tissue  density  in  grams  of 
tissue  per  cubic  centimeter  of  tissue,  and  the  specfic  he^t  of  the  blood 
(typically  b = .0122), 

and 

ITIS(I)  =*  the  tissue  type  indicated  by  a positive  integer  where  ITIS(I)  = 
1,2, 3, 4, 5, 6,  or  7 denotes  cerebrospinal  fluid,  blood,  muscle,  skin  or  dura, 
brain,  fat  or  bone,  or  yellow  bone  marrow,  respectively. 
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If  EPSP{I)  or  SIGP(I)  are  read  In  as  O.DO,  then  the  values  of  EPSP(I)  or 
SIGP(I)  or  both  are  replaced  by  values  determined  from  values  stored  in  data 
tables  by  means  of  the  commands: 

55  IF(EPSP(I).NE.O.DO.AND.SIGP(I).NE.O.DO) 

IGO  TO  60 

CALL  EPROP(FREO,ITIS(I),EP,SIG) 

IF(EPSP(I).EQ.O.DO)  EPSP(I)  = EP 
IF(SIGP(I).EQ.0.D0)  SIGP(I)  = SIG 

60  FAC2  = EPSP(I)/2.D0 

The  fourth  data  set  contains  NOCR  cards  each  containing  the  spherical 
coordinates  (R,THETAD,PHID)  of  a point  whose  cartesian  coordinates  are  (X,Y,Z) 
= (R*SIN(THETAD)*COS(PHID),R*SIN(THETAD)*SIN(PHID).R*COS(PHID))  at  which  the 
microwave-induced  temperature  rise  is  to  be  computed.  The  cards  are  read  in 
via  the  statements: 

30  F0RMAT{7F10.0) 

DO  345  1 1=1, NOCR 
READ  30,R,THETAD,PHID 
345  CONTINUE 

Finally  a typical  problem  is  presented  as  it  would  be  given  to  the  user 
and  tne  proper  response  is  indicated.  The  user  directions  are  to  compute  the 
thermal  response  of  a one-layer  spherically  symmetric  ball  of  brain  tissue 
with  a 2,804-cm  radius,  a permittivity  of  31,09,  a microwave  conductivity  of 
,0012,  a density  of  1,0,  a specific  heat  of  ,84,  and  a blood  flow  perviousness 
term  of  0,0  to  a steady  30-s  exposure  to  2450-HHz  radiation  with  a power  of  70 
mW/cm^.  The  user  is  to  compute  the  temperature  along  the  x,  y,  and  z-axes 
with  X,  y,  and  z-values  being  taken  ^rom  the  collectior:,  -2,8,  -2,45,  -2,10, 


-1.75,  -1.40,  -1.15,  -.8,  -.45,  -.1,  -l.E-4,  +1.E-4,  +.1,  +.45,  +.8,  +1.15, 
+1.40,  +1.75,  +2.10,  +2.45,  and  +2.8.  The  user  is  to  obtain  these  results 
with  maximum  accuracy.  The  proper  response  is  indicated  in  Figures  4.4.2  - 

4.4.5. 

2804.0-3  3109.D-2  1414.0-3  12.D-4  l.DO  .840+0  O.DO 

30.00  30.00  30.00  30.00  30.00  1 60 

2450.00  70.00  1.035  1 12  25  80  48  1 0 


Figure  4.4.2.  The  first  three  data  sets  for  the  computation  of  the  thermal 
response  of  a one-layer  brain  tissue  structure  exposed  to  70 
mW/cm^  and  2450-MHz  radiation  for  30  s at  60  si  points. 
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Figure  4.4.3.  Oata  set  describing  points  on  the  z-axis  in  spherical 
coordinates.  The  points  on  the  z-axis  at  which  the 
temperature  is  to  be  computed  are  shown.  The  columns  in  which 
the  data  entries  end  are  respectively  10,  20,  and  30. 
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Figure  4.4.4. 

Data  set  describing  points  on  the  x-axis  in  spherical 
coordinates.  The  points  on  the  x-axis  at  which  the  data 
entries  end  are  shown.  The  columns  in  which  the  data  entries 
end  are  respectively  10,  20,  and  30. 
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Figure  4.4.5. 

Data  set  describing  points  on  the  y-axis  in  spherical  coordi- 
nates. The  points  on  the  y-axis  at  which  the  temperature  is  to 

be  computed  are  shown.  The  columns  in  which  the  data  entries 
end  are  respectively  10,  20,  and  30. 
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4.5.  The  Output  and  its  Meaning 

The  output  of  our  program  to  predict  the  thermal  response  of  a spheri- 
cally symmetric  body  to  microwave  radiation  includes  (i)  the  printing  of  the 
input  data  defining  the  scattering  problem,  (ii)  the  weight  and  volume  of  the 
scatterer,  (iii)  the  average  and  total  absorbed  power,  (iv)  the  eigenvalues 
associated  with  radial  eigenfunctions,  (v)  the  expansion  coefficients  used  in 
expanding  the  temperature  in  spherical  harmonics,  and  finally  (vi ) the  pre- 
dicted microwave-induced  temperature  Increases  and  estimates  of  the  theoreti- 
cal error  in  our  predictions. 

The  input  data  which  defines  the  input  radiation  is  printed  out  and  is 
identified  by  labels.  This  data  includes  the  frequency,  the  field  strength, 
STOPR  which  is  defined  in  Section  4.3,  the  up  time  of  a single  pulse  or  TDUR, 
its  period  or  TPER,  the  number  NPUL  of  pulses  in  a single  pulse  train,  the 
time  TBPER  that  a single  pulse  train  lasts  which  includes  the  quiet  period 
after  the  initial  burst  of  NPUL  pulses,  the  time  TOUT  after  which  the  incident 
wave  is  cut  off  and  the  TIME  at  which  one  observes  the  temperature;  this 
output  data  set  is  printed  by  the  commands: 


PRINT  20,  FREQ,E01,UNIT(IE0+1),ST0PR,TDUR,TPER, 
1NPUL,TBPER,TCUT,TIME 

20  FORMATCITHERMAL  RESPONSE  OF  CONCENTRIC  SPHERICAL', 
I'HEAD  MODEL  TO  RFR '/'-FREQUENCY  =',F9.2, 

1'MHZ FIELD  STRENGTH  =' ,F9.2,1X,A8, 

16X, 'STOPR  ='1PD12.4/ 

I'O  FOR  ONE  PULSE,  UP  TIME  IS' ,012.4, 'SEC 

1'  AND  PERIOD  IS',D12.4,'SEC.7 

I'O  ONE  PULSE  TRAIN  CONTAINS' ,14, 'PULSES  AND  LASTS', 

lD12.4,'SEC.'/’0  THE  INCIDENT  WAVE  IS  CUT  OFF  AFTER', 

1D12.4,'SEC.  AND  THE  TEMPERATURE', 

1 'IS  OBSERVED  AFTER' ,D12.4, 'SEC. ' ) 
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The  next  set  of  output  data  defines  the  body  in  which  the  microwave- 
induced  temperature  is  to  be  predicted.  We  print  out  NORG  lines  of  data,  each 
of  which  defines  the  outermost  bounding  sphere  of  radius  SBDP(I),  the  relative 
permittivity  EPSP(I),  the  microwave  conductivity  SIGP(I),  the  thermal  conduc- 
tivity TCP(I),  the  density  RH0P{I),  the  specific  heat  CP{I),  the  blood  flow 
radial  perviousness  term  BFRP(I),  and  the  tissue  type  ITIS(I)  for  the  Ith 
tissue  layer.  This  output  data  set  is  defined  by  the  following  lines: 


PRINT  25 

25  FORMAT ( '-REGION* ,3X, 'SURFACE ' ,4X, 'RELATIVE' ,5X, 
l'ELECTRIC',7X. 'THERMAL', 6X, 'DENSITY', 3X, 

1 ' SPEC I F I C ' , 3X , ' BLOOD  FLOW ' , 4X , ' T I SSUE ' /9X , 

1 ' BOUNDARY ' , 3X , ' DIELECTRIC ' , 2X , ' CONDUCTIVITY ' , 

12X , ' CONDUCT I V ITY ' , 16X . ' HEAT ' , 8X , ' RATE ' . 7X , 

1 'TYPE' /21X, 'CONSTANT '/IIX, ' (CM) ' ,20X, ' (MHO/M) ' ,3X, 
1 ' (CAL/CM-SEC-C) ' .3X, ' (G/CM3) ' ,2X, ' (CAL/G/S) ' ,4X, 
l'(CC/SEC)'/) 


00  65  I = l.NORG 

65  PRINT  70,I,SBDP(I),EPSP(I),SIGP(I),TCP(I).RH0P(I), 
1CP{I),BFRP(I).TISSUE(ITIS(I)) 

70  F0RMAT(I4.F12.2,F13.3,F15.6,F13.4,F10.3, 

1F12.5,3X,A8) 

The  next  set  of  output  data  describes  intermediate  ouLpul  resulting  from 
defining  the  microwave  heat  source  term  for  the  heat  transfer  equation.  The 
data  printed  out  includes  the  mass  XMASS,  in  kilograms,  of  the  scattering 
body,  its  volume  VOL,  in  cubic  meters,  the  average  absorbed  power  PAVG  per 
unit  volume  expressed  in  watts  per  cubic  meter,  and  the  average  absorbed  power 
PAVGl  per  unit  mass  expressed  in  milliwatts  per  kilogram.  This  mode  by  which 
output  is  printed  is  described  by  the  statements: 
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PR I NT  80 , XMASS , VOL , PAV6 , ZLAB ( 2 ) , PAVGl , ZLAB ( 1 ) , TOTPOW 
80  F0RMAT('-WEIGHT=MPD12.4,'KG7'0  VOLUME*' ,D12. 4, 
l'M**3'/'0F0R  A CONTINUOUS  WAVE  THE  AVERAGE', 

1 'ABSORBED  POWER  IS' ,D13.5,A7, 'OR' ,D13.5.A7/ 

I'OTOTAL  ABSORBED  POWER*' ,D13. 5, 'WATTS'/ 

1' -ROOTS  OF  THE  EIGENFUNCTION'/) 

The  last  phrase  '-ROOTS  OF  THE  EIGENFUNCTION'  Is  a heading  for  the  next  output 
data  set  which  is  the  set  of  eigenvalues  needed  to  define  the  eigenfunctions 
used  in  expressing  the  microwave-induced  temperature  excursion. 

The  next  set  of  output  data  provide  us  with  the  array  XLAMDA  of  eigen- 
values defined  by  equations  (3.3.42)  and  (3.3.44),  and  which  are  used  to 
define  the  radial  eigenfunctions  used  in  expanding  the  microwave-induced 
temperature  excursion.  The  eigenvalues  are  printed  using  the  statements: 

DO  90  NSBF  * 1,NMAX 
Nl”  NSBF  - 1 

90  PRINT  95,N1,(XLAMDA(K,NSBF),K  * l.kMAX) 

95  F0RMAT(1H  ,I5,1PD12.4,9PD12.4/(7X,10D12.4)) 

Each  row  of  printing  in  this  output  data  set  displays  the  sequence  defined  by 
equation  (3.3.43)  where  the  row  index  N1  denotes  the  actual  Bessel  function 
order  and  K is  the  index  of  the  sequence  in  equation  (3.3.43). 

Once  the  eigenvalues  defined  as  the  solution  of  equation  (3.3.44)  are 
determined,  we  can  compute  the  radial  transform,  defined  by  equation  (3.3.22), 
of  the  source  term  and  subsequently  obtain  the  expansion  coefficients 
U(NSBF,M2,K),  defined  by  equation  (3.4.2),  that  are  used  in  representing  the 
microwave-induced  temperature  TRM.  The  expansion  coefficients  are  indexed  by 
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NSBF , which  is  one  plus  the  order  of  the  Bessel  function  under  consideration, 
M2  which  is  1 if  the  order  of  the  cosine  transform  is  0 and  is  2 if  this  is 
the  order  of  the  associated  cosine  transform,  and  finally  K which  is  the  index 
of  the  eigenvalue  associated  with  a given  Bessel  function  order.  The  expan- 
sion coefficients  are  printed  out  through  the  instructions: 


DO  270  NSBF  = l.NMAX 


DO *270  M = l.NSBF 


00*260  NRT  = l.KMAX 


260  LI(NSBF.M2.NRT)  = ETIME(NRT)*F*XNUM(NSBF,M2,NRT)/ 
1DEN(NSBF,M2,NRT) 

PRINT  265,N1.M1,(U(NSBF,M2,K),K  * l.KMAX) 

265  F0RMAT(2I3,1P10D12.4/(8X,1CD12.4)) 

270  CONTINUE 


The  final  and  most  important  output  describes  the  location  of  the  point  at 
which  the  temperature  is  sought,  the  predicted  microwave-induced  power  den- 
sity, the  temperature  excursion,  and  an  estimate  of  the  error  associated  with 
approximation  of  the  infinite  sum,  defined  by  equation  (3.4.1),  by  only  a 
finite  sum.  This  output  is  described  by  the  statements: 


PRINT  275,ZLAB{ISAR+1) 

275  F0RMAT{'r',29X' INTERNAL  POINT* ,11X.' ABSORBED ' , 
1 ' POWER ' ,7X, 'TEMPERATURE ' ,8X, 'APPROXIMATE ' , 
l/llX/'POINT' .2X, 'REGION' ,2X, 'RADIUS' ,3X. 

1 'THETA' ,4X, 'PHI ' ,12X, 'DENSITY' ,14X, 

1 ' R I SE • . 14X , 'ERROR ' /28X , ' CM ' ,6X , ' DEG ‘ , 12X , 
1A7,13X,’0EG  C',13X,'PER  CENT'/) 


DO  345  II  = l.NOCR 


PRINT  340, T r;E6,SAVR,THETA0,PHID,PD, 

1TRM,PCEBF,PL.., 

340  F0RM/a(I14,I8,F10.3,2F8.2,F19.8,lPD20.4,2P2F12.: , 

345  CONTINUE 

In  tha  above  2LAB(ISA.»+-1)  is  an  alpha  array  cc’taining  the  label  'MW/KG'  which 
;■^3r^ds  for  milliwatts  per  kilogram  or  'W/M**3'  which  stands  for  watts  per 
'•'^c  mete''.  The  parameters  II  and  NREG  denote,  respectively,  the  index, 

r.  1 ng  .'rotr;  > r NOCR,  of  the  poinc  at  which  the  temperature  is  to  be 
cci..p''ted  an'”  t.  >:  layer  number,  ranging  from  1 to  NORG,  of  the  layer  in  which 
the  point  is  found.  The  three-tuple  (SAVER,  THETAD,  PHIO)  is  the  spherical 
coord'-  ’•e  representation  of  the  point  at  which  the  temperature  is  to  be 
computed.  T*!?  variables  PD  and  TRM  denote,  respectively,  the  microwave  power 
per  unit  volume,  and  the  microwave-induced  temperature  excursion  at  the  point 
(SAVER,  THETAD,  PHID),  The  error  estimation  parameter  PCEBF  denotes  the 
relative  error  in  temperature  prediction  associated  with  using,  NMAX  - 1 
orders  of  Btssel  functions  and  KMAX  eigenvalues  per  Bessel  function  instead  of 
using  NMAX  and  KMAX  to  get  a temperature  estimate  SBFMl.  On  the  other  hand, 
we  see  whether  or  not  we  have  used  enough  eigenvalues  per  Sessel  function 
order  by  using  NMAX  orders  of  Bessel  furctioris  and  KMAX  - 1 eigenvalues  per 
Bessel  functions  order  obtaining  : temperature  estimate  SRMl  and  computing  the 
relative  error  PCER  by  the  statement, 

PCER  = (TRM  - SRMl) /TRM 


4,6.  Program  Size  and  Running  Time 

The  program  requires  252K  on  the  'GO'  step  for  an  IBM  360  and  has  a run- 
ning time  that  is  dependent  on  the  accuracy  demanded  and  the  number  of  layers 
in  the  model.  For  a one-layer  model  demanding  maximum  accuracy  and  computing 
the  temperature  at  60  points  for  one  exposure  time,  the  time  on  the  'GO'  step 
was  2.93  minutes. 

Gaussian  quadrature  i:  used  to  compute  cosine,  Legendre,  and  radial 
transforms  of  the  source  ta>ni  divided  by  the  product  of  the  density  and 
specific  heat.  We  do  these  compulations  in  an  optimal  way  by  precomputing 
needed  values  and  taking  care  not  to  compute  any  complex  function  more  than 
once  at  the  same  argument. 


4.7,  Error  Messages 


In  this  section  we  explain  the  error  messages  that  the  program  provides 
to  the  user  when  he  has  inadvertently  provided  unsuitable  input  data.  Some  of 
the  errors  are  fatal  and  some  merely  provide  a warning  to  the  user  regarding 
the  accuracy  of  their  results. 

An  example  of  the  latter  occurs  often  when  one  attempts  to  compute  the 
thermal  response  at  points  on  the  positive  z-axis  in  that  the  series  expansion 
of  the  electric  field  vector  may  not  have  enough  terms  in  it  to  guarantee 
eight  significant  digits  of  accuracy  in  the  answer.  The  coding  which  prints 
out  this  error  message  is  given  by  the  statements: 

PRINT  30,NMIN,NC,ST0PR,EPS 

30  FORMAT (15X,'NMIN  ='.13,'  NC  ='.I3.2X. 

I'STOPR  =',1PD14.4,'IS  TOO  SMALL', 

I'FOR  ACCURACY  0F',D14.4) 

where  NC  is  the  number  of  spherical  Bessel  functions  available  to  estimate  the 
field  at  a given  point,  NMIN  is  the  number  of  expansion  coefficients  available 
based  on  the  location  of  the  layer  electrical  properties,  and  EPS  is  the 
relative  error  demanded  in  the  solution. 

The  error  messages  are  described  next  in  the  order  in  which  they  are 
found  in  the  main  program.  The  first  message  in  the  main  program  gives  an 
obvious  constraint  on  the  parameters  defining  the  time  profile  of  the  beam. 
This  is  printed  out  when  appropriate  by  the  statements: 

IF  (NPUL.GT.O.AND.TDUR.GT.O.DO. 

1AND.TCUT.GT.0.D0.AND.TIME.GT.0.D0)G0  TO  24 

21  PRINT  22 

22  FORMAT f'****ERR0R  IN  TIMES****') 
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24  IF(TPER.LT. TOUR.OR.TBPER.lt. NPUL*TPER)GO  TO  21 


The  next  fatal  error  messages  deal  with  the  fact  that  the  radii  of  the 

layers  should  be  in  ascending  order  and  that  there  are  only  7 possible  tissue 

types,  recognized  by  the  program,  assignable  to  a layer.  These  messages,  when 

appropriate,  are  printed  by  the  commands: 

IF(I.EQ.1.0R.SBDP{I).GT.SBDP(I-1))60  TO  45 
PRINT  40 

40  FORMAT ('****L AYER  RADII  MUST  BE'. 

1'  IN  ASCENDING  ORDER****) 

STOP 

45  IF(ITIS{I).GT.0.AND.ITIS(I).LT.8)60  TO  55 
PRINT  50,I,ITIS(I) 


The  next  control  on  the  input  is  based  on  the  fact  that  the  number  of 

points  used  in  the  Gauss  guadrature  Integration  scheme  for  determining  the 

expansion  coefficients  can  only  assume  certain  discrete  values  contained  in  a 

five-element  array  NPOINT.  These  error  messages  controlling  the  number  of 

points  requested  to  be  used  in  evaluating  the  radial  transform  are,  when 

appropriate,  printed  by  means  of  the  commands: 

DO  100  I = 1,5 
IF(MP.EQ.NPOINT(I))GO  TO  110 
100  CONTINUE 
PRINT  105,MP 

105  FORMATCOINTEGRATION  CONTROL  =',I9,2X, 

I'IS  NOT  AVAILABLE') 


The  error  message  controlling  the  number  of  points  requested  to  he  used 
in  evaluating  the  Legendre  transform  is,  when  appropriate,  printed  by  means  of 
the  commands: 


DO  115  I = 1,5 

IF{MP1.EQ.NP0INT(I))G0  TO  120 
115  CONTINUE 

PRINT  105, MPl 
STOP 


The  above  errors  are  fatal  in  the  sense  that  the  program  stops  execution 
as  soon  as  the  errors  are  recognized. 

The  next  error  message  is  another  check  on  the  accuracy  with  which  the 
electric  vectors  are  computed.  This  check  deals  with  the  computation  of 
fields  at  the  Gauss  quadrature  points  for  the  purpose  of  eliriuating  the  tem- 
perature excursion  induced  by  the  microwave  radiation.  The  message  described 
at  the  beginning  of  this  section  dealt  with  the  computation  of  power  density 
at  the  points  at  which  the  temperature  is  to  be  computed  or,  said  differently, 
with  the  accuracy  of  column  6 in  the  last  output  data  table. 

These  error  messages  are  described  by  the  following  statements: 

IR  = 0 

DO  165  N = 1,NC 
FACl  = 2*N  + 1 
IF(IR.E0.1)GO  TO  155 
T = P(N)*TR(N) 

ERAD  = ERAD  + T 

IF(CDABS(T).LT.CDABS(ERAD)*EPS)  IR  = 1 
155  IF{ITP.E0.1)GO  TO  160 
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NPl  = N + 1 
RATIO  = FAC1/(N*NP1) 

A = RATIQ*P(N)/SINTH 
B = -RATIO*DP(N) 

C = A*TEI(K)  + B*TE(N) 

ETHETA  = ETHETA  + C 
T = A*TEI(N)  + B*TE(N) 

EPHI  = EPHI  + T 

1F(CABS(C).IT.CDABS(ETHETA)*E?S. 
1AND.CDABS{T).LT.CDABS(EPHI)*EPS)  ITP  = i 
160  IF{IR  + ITP.EQ.2)G0  TO  175 
165  CONTINUE 

PRINT  1 70 , NMI N, NC .THETA ,R .STOPR .EPS 
170  F0RMAT(15X.'NMIN='I3.'NC='.I3. 
l'THETA='.F9.6.'R='.2PF9.6. 
1'ST0PR='1P09.2.'IST00  SMALL  FOR*. IX. 
I'ACCURACY  0F'.D9.2) 

175  ERAD  = ERAD/Q 


We  note  that  in  the  above  code  C represents  an  amount  to  be  added  to  the 
series  representation  of  ETHETA,  Thus,  if  CDABS(C)  is  small  in  comparison  to 
CDABS(ETHETA)*EPS,  then  we  say  that  adding  the  term  C affected  the  value  of 
ETHETA  in  the  decimal  place  equal  to  the  integer  value  of  1/EPS  or.  said 
differently,  that  EPS  is  the  relative  error  associated  with  using  one  less 
term  in  the  series  representation  of  ETHETA.  It  is  also  clear  from  the  above 
coding  that  the  term  T is  used  in  the  same  w^^y  to  describe  the  accuracy  with 
which  ERAD.  the  component  of  the  electric  field  vector  in  the  radial  direc- 
tion. and  EPHI.  the  component  of  the  electric  field  vector  in  phi  direction, 
are  computed. 

The  final  error  message  of  the  control  program  is  a nonfatal  error 
message  that  warns  the  user  when  he  attempts  to  predict  the  microwave  heating 
in  the  free  space  outside  the  body  that  is  being  irradiated.  This  message, 
when  appropriate,  is  printed  by  means  of  the  commands: 


106 


DO  285  NREG  = l.NORG 
IF(R.LE.SBDP(NREG))60  TO  300 
285  CONTINUE 
290  NREG  = 1000000000 

PRINT  295,U,NREG,SAVR,THETA0,PHID 
295  FORMAT{I14,I8.F10.3,3X,'**THE  RADIUS'. IX, 

I'IS  OUTSIDE  THE  SPHERE***') 

. GO  TO  345 
300  I REG  = NREG 

The  statement 

GO  TO  345 

directs  the  program  around  the  temperature  computation  part  of  the  orogram 
when  this  error  message  is  printed.  In  other  words,  the  computer  program  will 
not  let  the  user  waste  his  time  by  attempting  to  compute  microwave-induced 
temperature  excursions  at  points  outside  the  body  being  irradiated. 


4.8.  Program  and  Subprogram  Oescription 

In  this  section  we  give  an  executive  description  of  the  overall  program, 
list  the  subroutines  called,  and  give  their  purpose. 

The  main  program  is  divided  into  five  parts.  These  parts  carry  out  (i) 
the  scattering  problem  definition  by  reading  in  the  data  and  using  subroutine 
EPROP,  (ii)  the  electromagnetic  field  expansion  coefficient  determination  and 
surface  energy  balance  from  the  results  of  the  COEF  subroutine,  (iii)  the 
determination  of  the  eigenvalues  of  the  elliptic  part  of  the  heat  transfer 
operator  using  the  RFNDR  subroutine,  (iv)  the  microwave  heat  source  expansion 
and  thermal  expansion  coefficient  using  the  subroutines  BJYH,  TERM,  PL,  ALP, 
and  SRBF,  and  finally  (v)  the  power  density,  temperature  and  error  estimation 
portion  using  only  the  subroutines  SRBF  and  ALP.  The  beginning  and  ending  of 
the  above  five  sections  are  marked  by  comment  cards  in  the  listing  of  the 
program  in  Appendix  A. 

In  the  next  part  of  this  section  we  describe  all  of  these  subroutines  in 
the  order  in  which  they  occur  in  the  main  f,rogram. 

The  subroutine  EPROP  determines,  by  interpolating  tabulated  data,  the 
relative  permittivity  or  microwave  conductivity  of  any  of  the  seven  tissue 
types  from  tabulated  data.  The  subroutine  is  called  by  the  statement, 

CALL  EPR0P(FREQ,ITIS{I),EP,Si8) 

The  user  must  supply  FREQ,  the  frequency  of  the  incoming  radiation  in  mega- 
hertz, and  the  tissue  type  ITIS{I)  of  the  Ith  layer  of  the  scatterer,  where 
ITIS(I)  = 1,  2,  3,  4,  5,  6,  or  7 depending  on  whether  the  tissue  type  is 
cerebrospinal  fluid,  blood,  muscle,  skin  or  dura,  brsin,  fat  or  bone,  or 
yellow  bone  marrow,  respectively. 
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The  subroutine  COEF  generates  expansion  coefficients  ANP,  BNP,  ALPNP,  and 
BETNP  used  in  expanding  the  electromagnetic  field.  It  is  called  by  the  state- 
ment, 

CALL  COEF 

The  subroutine  RFNDR  is  used  to  determine  the  eigenvalues  of  the  elliptic 
part  of  the  heat  transfer  operator  by  a shooting  method.  Basically  we  start 
out  with  a trial  value  of  the  eigenvalue,  an  assumption  about  the  asymptotic 
behavior  of  the  radial  eigenvalue  at  the  origin  so  that  with  this  assumption 
the  solution  of  the  singular  ordinary  differential  equation  with  which  the 
eigenvalue  is  associated  is  unique.  Wa  then  check  and  see  if  the  Newton  cool- 
ing condition  on  the  boundary  is  satisfied.  If  it  is,  we  know  that  we  have  an 
eigenvalue.  If  the  Newton  cooling  condition  or  equation  (3.3.44)  is  not 
satisfied,  we  increase  the  trial  value  slightly  and  try  again. 

When  we  find  two  trial  values  at  which  the  Newton  cooling  function,  the 
output  of  the  function  subroutine  FNCAL,  differ  in  sign,  we  then  use  a bisec- 
tion routine  DRTMI  to  get  the  value  of  the  eigenvalue  to  as  many  decimal 
places  as  is  desired. 

The  subroutine  BJYH  generates  arrays  BJNP  and  BHNP  of  spherical  Bessel 
functions  and  Hankel  functions,  respectively,  used  In  the  determination  of  the 
functions  defined  in  equations  (3.2,2)  - (3.2.6).  It  is  called  by  the  state- 
ment, 

CAL  BJYH(BJNP,BHNP,Q,NC,STOPR,HAX) 

where  0 is  equal  to  a sphere  radius  multiplied  by  the  complex  propagation 
constant  FKP(NREG)  defined  by  equation  (3.2.7).  We  compute  up  to  MAX  values 
limited  by  the  size  constraint  STOPR,  but  we  fill  the  arrays  with  only  NC 
values  since  we  have  the  same  number  of  spherical  Bessel  functions  in  all 
regions. 
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The  subroutine  TERM  computes  the  product  of  the  square  root  of  (-1) 
raised  to  the  power  NCK  and  the  factor  of  I**L  appearing  in  equation  (3.2.2) 
where  here  "I"  denotes  the  square  root  of  -1. 

The  subroutine  PL  computes  an  array  P of  associated  Legendre  polynomials 
of  the  first  kind  and  order  1 and  an  array  DP  of  their  derivatives. 

The  function  subroutine  ALP  computes  the  associated  Legendre  function  of 
the  first  kind,  degree  N and  order  M with  M and  N being  nonnegative  integers 
and  with  M not  exceeding  N.  It  is  a function  subroutine  returning  a single 
value  at  a single  point. 

The  function  subroutine  SRBF  computes  the  spherical  Bessel  functions  XJ 
and  XY  of  the  first  and  second  kind,  respectively,  and  their  respective  deriv- 
atives DJ  and  DY  at  the  value  equal  to  the  square  root  of  SI  multiplied  by 
SAVR  or  the  arguments  appearing  in  the  discussion  in  Section  3.3. 
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PROGRAM  TRP  1 

THERMAL  RESPONSE  OF  CONCENTRIC  SPHERICAL  2 

HEAD  MODEL  TO  RFR  3 

4 

IMPLICIT  REAL*8  (A-H,0-Z)  5 

C0MPLEX*16  FKP,ANP,BNP,ALPNP.BETNP,BJNP,BHNP,ERAD,ETHETA,EPHI,T,C,  6 

1W,X,Y1,Z,0,TE{50),TE1(50),TR(50)  7 

COMMON  FKP(7),ANP(300),BNP(300).ALPNP(300).BETNP(300),BJNP(100),BH  8 

1NP(100),BDP(6),P(51),DP(50),R.THETA,C0STH,PHI,SINTH,ST0PR,E0  9 

COMMON  /A/NORG,NREG,NRT,NSBF.NMIN,NC,ICODE  10 

COMMON  /B/FACT(6,25,18),AJ(6,25.18),BY(6,25,18),XLAMDA(25,18),SBDP  11 

1(6),RH0P(6),CP{6).BP(6),TCP(6),H  12 

COMMON  /C/AJ1,S1,F,R1,IREG  13 

INTEGER*2  IFL(102,102)  14 

REAL *4  R3(304),TR3(304),X2(102),DAR(102,102),CLAB(3,3),ANG,AX1,AY  15 

DIMENSION  U{18,  2,25),EPSP(6).SI6P(6),BFRP(6),  16 

1 sr80,64,2),XNUM(18,2,25),DEN(18.2,25),RR(80),ETIME{25)  17 

1 . THET1{64),C0STH1(  18 

164),SINTH1(64),WTTH{64).ALP0L(64)  , 19 

1 NP0INT(5),KEY(6).Y(116),  WT(116),  ARG3(2)  20 

1,UNIT{2),ZLAB(2),  21 

1 SUM2(80,2)  , 22 

1 TISSUE(7),ITIS(6)  23 

1 ,BLAB(3),AX(3),  DLAB(4)  24 

DATA  TISSUE/'CS  FLUID' , 'BLOOD ‘ , ’MUSCLE' . 'SKIN-DUR 'BRAIN' , 'FAT-BO  25 
INE'.'Y.B.M.'/  26 

DATA  UNIT/'V/  27 

IM'  ,'MW/CM**2'/,ZLAB/'  MW/KG', ' W/M**3'/  28 

1 , EPS/1. 0-8/  29 

DATA  CLAB/'E  PL'.'ANE','  ', 'H  PL ',  ANE', ' ','X-Y','PLAN', 'E'/  30 

1,BLAB/'Z-AXIS  C','00RDINAT'  31 

1 ,'E  (CM)'/,AX/'Z-AXIS  C, 'X-AXIS  C. 'Y-AXIS  C'/,  32 

1 OLAB/'TEMPERAT','URE  RISE','  (D  33 

lEG.  C ' , ' . ) '/  34 

DATA  NPoiNT/32, 48,64, 80, 8/, KEY/1, 17,41, 73, 113, 117/  35 


DATA  Y/.048307665688D0, . 1444719615800, .2392873622500, .33186860228D  36 

10..  42135127613D0,.50689990893D0,.58771575724D0, .6630442669300,. 732  37 

11821187400..  7944837959700. .84936761373D0,.89632115577D0,. 934906075  38 

194D0,.96476225559D0,.98561151155D0,.99726386185D0,.0323801/0963D0,  39 

l.O97O04699209DO,.161222356O7DO,.22476379O39D0,.28736248736D0,.3487  40 

1 5588629D0 , .40868648199DO , .46690290475D0 , . 52316097472D0 , . 5772247260  41 
18D0, .6288673967800, .67787237963D0, .7240341309200, .7671590325200, .8  42 

1070662040300 . . 8435882616200 . .8765720202700 . .9058791367200 . . 9313866  43 

1907100..  9529877031600.. 9705915925500.. 9841245837200.. 9935301722700  44 

1 . .  99877100725D0 , .02435029266300 , .07299312178800 , . 1214628193D0 , . 169  45 

16444204200..  2174236437400.. 2646871622100.. 3113228719900.. 357220158  46 

13400..  4022701579600.. 4463660172500. .4894031457100. .5312794640200..  47 

1571895646200..  6111553551700.. 6489654712500.. 6852363130500.. 7198818  48 

1501700..  7528199072600.. 7839723589400.. 8132653151200.. 8406292962500  49 

1. . 8659993981500. .88931544600. .9105221370800. .9295691721300. .946411  50 

13748600 . .  9610087996500 . . 9733268277900 . . 9833362538800 . . 991013371480  51 

10..  9963401 167700.. 9993050417400,  52 

1 .01951138325700, .05850443715200, .09740839844200, .136164022  53 


I 


1 


I 

I 


j 
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1809D0 , . 1 7471229183D0 , .2129945028600,. 25095235839D0 , . 28852805488D0 , 54 

1.32566437075D0,. 362304753500,. 3983934058800,. 43387537083D0,. 468696  55 

16151700..  5028041118900.. 536145920900.. 5686712681200.. 6003306228300  56 

1..  6310757730500.. 6608598989900.. 6896376443400.. 7173651853600.. 7440  57 

1002975800 . .  7695024201400 . . 793832717500. .8169541386800 . .83883147358  58 

100. . 8594314066600. .8787225676800. .8966755794400. .9132631025700. .92  59 

184598771700 . .  9422427613100 . .9545907663400 . . 9654850890400 . . 97490914  60 

105900..  9828485727400.. 989291302500. .9942275409700.. 997649864400.. 9 61 

1995538226500 . .  183434642500 . . 5255324099200 . . 7966664774100 . . 96028985  62 

16500/  63 

DATA  WT/. 09654008851500,. 09563872007900,. 09384439908100,. 091173878  64 

169600. . 08765209300400. .08331192422700. .07819389578700. .07234579410  65 

1900 . . 06582222277600 . . 058684093479D0 , .05099805926200 , . 0428358980220  66 

10..  03427386291300.. 02539206530900.. 01627439473100.. 007018610009500  67 

1..  06473769681300.. 06446616443600. .06392423858500.. 06311419228600..  68 

10620394231600..  06070443916600.. 05911483969800.. 057277292100.. 05519  69 

19503700..  05239018948500.. 05035903555400. .04761665849200.. 044674560  70 

185700..  04154508294300.. 0382^135106600,. 03477722256500,. 03116722783  71 

1300..  02742650970800.. 02357076083900. .01961616045700.. 0155793157230  72 

10..  01147723457900.. 007327553901300. .003153346052300.. 0486909570090  73 

10 . .  04857546744200 . . 04834476223500 . .04799938859600 . . 04754016571500 , 74 

1 .04696818281600. . 04628479658100. .04549162792700. .04459055816400. .0  75 

14358372452900..  04247351512400. .04126256324300.. 03995574113300. .038  76 

155015317900. . 0370551285400. .03547221325700. .03380516183700. .032057  77 

192835500..  03023465707200.. 02833967261400.. 02637746971500. .02435270  78 

1256900..  02227017380800.. 02013482315400.. 017951715776D0,. 0157260304  79 

1 7600 . . 01346304789700 ..011 1681394600 . .008846759826400 . . 006504457969  80 

100. . 004147033260600. .001783280721700,  81 

1 .U39017813656D0,. 03895839596300, .03883965105900, .03866175  82 

1977400. . 03842499300700. .03812971131400. .03777636436200. .0373654902  83 

13900..  03689‘'71463800,. 03637374990600, .03579439395300, .035160529045  84 

100..  03447312045200,  03373321498500,. 03294193939800,. 03210049867300  85 

1.. 03121017418Bnn,. 0302723^17600, .02928836958300,. 02825981605700, .0  86 

1271882275CG, .02607523576800, .02492253576400, .02373188286600, .02250  87 

1509024600 . .  02 12440261 16D0 , .01995061087800 , . 01862681420800 , . 0172746  88 

15205600. . 01589618358400. .01449350804100. .01306876159200. .011624114  89 

112100..  01016176604100.. 00C6339452693D0,. 007192904768100,. 005690922  90 

1451400. . 004180313124700. .002663533589500. .001144950003200. .3626837  91 

1833800 . .  3 137066458800 . . 22238103445D0 , . 1012285362900/  92 

IPLSW=0  93 

PIE=3. 14159265358979300  94 

RAD=180. DO/PIE  95 

EPS0=8. 854160-12  96 

VEL=2.997924562D8  97 

RH0BP=1.0600  98 

CBP=0.98D0  99 

H=5. 720-5  100 

ITMfc'=0  101 

**********FIRST  DATA  CARD  - CONTROL  PARAMETERS  102 

READ  (5,5)FREO,EO,STOPR,NORG,NMAX,KMAX,MP,MP1,IEO,ISAR  103 

5 FORMAT  (3010.0,715)  104 

FREQ  FREQIIPNCY  IN  MEGAHERTZ  105 

EO  STRE'-  H OF  INCIDENT  E-FIELD  106 

STOPR  CUTOh.  FOR  SBF  COMPUTATIONS  107 

NORG  NUMBER  OF  LAYERS  IN  SPHERE  108 

IS  DESIRED.  A CARD  WILL  BE  READ  FOR  EACH  POINT.  109 

NMAX  NUMBER  OF  ORDERS  OF  BESSEL  FUNCTIONS  USED.  MAX=12  110 
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KMAX  NUMBER  OF  ROOTS  FOR  EACH  ORDER.  MAX=J?5  111 

MR  NUMBER  OF  POINTS  FOR  INTEGRATION  FOR  RADIUS.  32,  48,  112 

64  OR  80.  {,  IS  AVAILABLE  FOR  TEST  RUNS)  113 

MPl  NUMBER  OF  POINTS  FOR  INTEGRATION  FOR  THETA.  32  OR  48  114 

lEO  INPUT  EO  UNITS  115 

0 - VOLTS/METER  116 

1 - MILLIWATTS/SQUARE  CENTIMETER  117 

ISAR  OUTPUT  POWER  DENSITY  UNITS  118 

0 - MILLIWATTS/KG  119 

1 - WATTS/CUBIC  METER  120 

• E01=E0  121 

IF  (JEO.EQ.O)  GO  TO  iO  122 

IE0=1  123 

E0-^DSQRT(3767.D0*E0)  124 

**********SECOND  DATA  CARD.  TIMES  IN  SECONDS  FOR  INCIDENT  WAVE  PULSES.  125 
FIRST  PULSE  TURNS  ON  AT  ZERO  SECONDS.  126 

10  READ  {5.15,END=495)  TDUR,TPER.TBPER.TCUT,TIME,NPUL,N0CR,IPL1 ,IPL2,  127 
INTR  128 

15  FORMAT  (5D10. 0,515)  129 

TDrR  TIME  DURATION  OF  A PULSE  130 

TPER  PERIOD  FROM  START  OF  A PULSE  TO  START  OF  NEXT  PULSE.  131 

TBPiIR  PERIOD  FOR  A GROUP  OF  PULSES.  132 

TCUT  TIME  AT  WHICH  WAVE  IS  CUT  OFF  133 

TIME  TIME  WHEN  TEMPERATURE  RISE  IS  OBSERVED.  134 

NPUL  NUMBER  OF  PULSES  IN  A GROUP  135 

NOCR  NUMBER  OF  POINTS  IN  SPHERE  AT  WHICH  TEMPERATURE  RISE  136 

***  PRINT  OUT  TITLE  AND  BASIC  INPUT  DATA  137 

PRINT  20,FR£Q,E01,UNIT(IE0+1),ST0PR,TDUR,TPER.NPUL,TBPER.TCUT,TIME  138 

20  FORMAT  ('ITHERMAL  RESPONSE  OF  CONCENTRIC  SPHERICAL  HEAD  MODEL  TO  R 139 

IFR '/'-FREQUENCY  =',F9.2,'  MHZ  FIELD  STRENGTH  «',F9.2,1X,A8, 6 140 

1X,'ST0PR  =',1PD12.4/  141 

rOFOR  ONE  PULSE,  UP  TIME  IS ',012.4,'  SEC.  AND  PERIOD  IS',D12.4,'  S 142 
2EC.7'OONE  PULSE  TRAIN  CONTAINS' ,14, ' PULSES  AND  LASTS' ,D12. 4, ' SE  143 
3C.  7'0THE  INCIDENT  WAVE  IS  CUT  OFF  AFTER' ,D12.4, ' SEC.  AND  THE  TEM  144 
4PERATURE  IS  OBSERVED  AFTER' ,D12.4, * SEC.')  145 

IF  (NPUL.GT.O.AND.TDUR.GT.O.DO.ANO.TCUT.GT.O.DO.AND.TIME.GT.O.DO)  146 
IGO  TO  24  147 

21  PRINT  22  148 

22  FORMAT  ('****  ERROR  IN  TIMES  ****')  149 

STOP  150 

24  IF  {TPER.LT.TDUR.OR.TBPER.LT.NPUL*TPER)  GO  TO  21  151 

ITME=ITME+1  152 

IF  (ITME.ST.l)  GO  TO  215  153 

PRINT  25  154 

25  FORMAT  ('OREGION  SURFACE  RELATIVE  ELECTRIC  THERMAL  155 

1 DENSITY  SPECIFIC  BLOOD  FLOW  TISSUE'/  156 

1 9X, 'BOUNDARY  DIELECTRIC  C 157 

lONDUCTTVITY  CONDUCTIVITY' ,16X, 'HEAT' ,8X, 'RATE  TYPE'/  158 

1 21 X,' CONSTANT '/I IX  159 

1,'{CM)',20X.'(HM0/M)  (CAL/CM-SEC-C)  {G/CH3)  (CAL/G/$)  (CC/  160 

ISEC)'/)  161 

0MEGA=2.D6*PIE*FREQ  162 

FAC1=0MEGA/VEL  163 

START=*1.D38  164 

XMASS=O.DO  165 

OLOVOL=O.DO  166 

**********REAO  LAYER  PROPERTIES  AND  COMPUTE  PROPAGATION  CONSTANTS  167 
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DO  65  I=l,NORG  168 

READ  (5,30)SBDP{r),EPSP(I),SIGP(I).TCP{I),RH0P(I),CP(I),BFRP(I),  169 
1 ITJS(I)  170 

30  FORMAT  (7F10.0,I5)  171 

SBDP  OUTER  RADIUS  OF  LAYER  IN  CENTIMETERS  172 

BDP  LAYER  OUTER  RADIUS  IN  METERS  173 

EPSP  PERMITTIVITY(RELATIVE)  174 

SIGP  CONDUCTIVITY  (MHOS  PER  METER)  175 

TCP  THERMAL  CONDUCTIVITY  176 

RHOP  DENSITY  177 

rP  '^PFf'TFir  HFAT  17fl 

BFRP  BLOOD  FLOW  RADIAL  PERVICIVITY  179 

ITIS  CODE  FOR  LAYER  TISSUE  TYPE  180 

1 DENOTES  CEREBROSPINAL  FLUID  181 

2 DENOTES  BLOOD  182 

3 DENOTES  MUSCLE  183 

4 DENOTES  SKIN  OR  DURA  184 

5 DENOTES  BRAIN  185 

6 DENOTES  FAT  OR  BONE  186 

7 DENOTES  YELLOW  BONE  MARROW  187 

BDP(I)=SBDP(I)/1.D2  188 

V0L=4.D0*PIE*B0P(I)**3/3.D0  189 

RV0L=V0L-0LDV0L  190 

XMASS=XMASS+RVOL*RHOP { I )*1 .03  191 

OLDVOL=VOL  192 

BP(I)=RHOBP*CBP*BFRP(n  193 

A=BP(I)/(RHOP(I)*CP(I))  194 

IF(A.LT.START)  START=A  195 

IF  (I.EQ.l.OR.SBDP(I).GT.SBDP(I-l))  60  TO  45  196 

PRINT  40  197 

40  FORMAT  ('0****  LAYER  RADII  MUST  BE  IN  ASCENDING  ORDER  ****•)  198 

STOP  199 

45  IF  (ITIS(I).GT.0.AND.:TIS(I).LT.8)  GO  TO  55  200 

PRINT  50,I,ITIS(I)  201 

50  FORMAT  ( •0****TISSUc  TYPE  CODE  FOR  LAYERM2,'  ISM5,’,  OUTSIDE  T 202 

IHE  RANGE  1-7****')  203 

STOP  204 

55  IF  (EPSP(I).NE.O.OO.AND.SIGP(I).NE.O.DO)  GO  TO  60  205 

CALL  EPROP(FREO,ITIS(I),EP,SIG)  206 

IF  (EPSP(I).EQ.O.DO)  EPSP(I):=EP  207 

IF  {SIGP(I).EQ.O.DO)  SIGP(I)=SIG  208 

60  FAC2=EPSP{I)/2.D0  209 

FAC3=DSORT{1.DO+(1.00/(EPSO*OMEGA)**2)*(SIGP(I)/EPSP(I))**2)  210 

FKP(I)=DCMPLX(FACl*DSQRT{FAC2*(FAC3+l.D0)),FACl*nSQRT(FAC2*(FAC3-l  211 
l.DO)))  212 

65  PRINT  70,I,SB0P(I),EPSP(I),SIGP(I).TCP(I),RH0P{I),CP(I).BFRP(I),TI  213 
1SSUE(ITIS{I))  214 

70  FORMAT  (I4,F12.2,F12.2,F13.3,F15.6,F13.4,F10.3,F12.5,3X.A8)  215 

FKP(N0RG-!l  )-0CMPLX  (FACl  ,0.D0 ) 216 

IF  (START. EO.O. DO)  START=l.D-9  217 

COMPUTE  EXPANSION  COEFFICIENTS  AND  TOTAL  ABSORBED  POWER  218 

CALL  COEF  219 

NN=NORG*NMIN  220 

QS^O.DO  221 

QT=O.DO  222 

DO  75  N=1,NMIN  223 

FACN=2*N+1  224 


117 


N?.!N=NN+N  225 

X3=FACN*nREAL (ALPNP (NNN }+BETNP (NNN ) ) 226 

Y3=FACN* (CDABS (ALPNP (NNN ) )**2+CDABS (BETNP (NNN ) )**2 ) 227 

QT=0T-X3  228 

QS=0S+Y3  229 

IF  (DABS(X3).LT.DABS(QT)*l.D-6.AND.Y3.LT.0S*l.D-6)  GO  TO  242  230 

75  CONTINUE  231 

PRINT  241  232 

241  F0RMAT('0****  TOO  FEW  EXPANSION  COEFFICIENTS  ****')  233 

242  TOTPOW=2.65441D-3*EO**2*2.DO*PIE*(QT-OS)/(2.DO*FAC1*FAC1)  234 

PAVG=TOTPOW/VOL  235 

PAVG1=1.D3*T0TP0W/XMASS  236 

***  PRINT  OUT  AVERAGE  ABSORBED  POWER  DENSITY  AND  TOTAL  ABSORBED  237 

***  POWER  238 

PRINT  80,XMASS,V0L,PAVG,ZLAB(2),PAVG1.ZLAB(1),T0TP0W  239 

80  FORMAT  ('OWEIGHT  =',1PD12.4,'  KG'/'OVOLUME  ='.D12.4,'  M**37  240 

rOFOR  A CONTINUOUS  WAVE  THE  AVERAGE  ABSORBED  POWER  IS' ,D13.5,A7. ' 241 

10R',D13.5,A7/  242 

I'OTOTAL  ABSORBED  POWER  =',D13.5.'  WATTS'/  243 

1 'OROOTS  OF  THE  EIGENFUNCTION'/)  244 

GET  ROOTS  OF  EQUATION  245 

STEP=l.D-8  246 

IC0DE=0  247 

AJ1=1.D0  248 

PRINT  543, START, STEP  249 

543  FORMAT  ('OSTART  =',1P015. 7,'  STEP  =' ,D15.7)  250 

DO  90  NSBF=1,NMAX  251 

AJ1=AJ1*(2*NSBF-1)  252 

CALL  RFNDR(START,STEP, 1.0-6, KMAX, 1000,10000)  253 

START=XLAMDA(1,NSBF)  254 

STEP= ( XLAMDA (2 , NSBF )-XLAMDA ( 1 ,NSBF ) )/10 .DO  255 

N1=NSBF-1  256 

PRINT  95, N1,(XLAMDA(K, NSBF), K«l, KMAX)  257 

95  F0RMAT(1H  ,I5,1PD12.4,9D12.4/(7X,10D12.4))  258 

IF  (NS8F.EQ.1)  GO  TO  90  259 

DO  89  K=1,KMAX  260 

IF  (XLAMDA(K,NSBF).LE.XLAMDA(K,NSBF-1))G0  TO  500  261 

89  CONTINUE  262 

90  CONTINUE  263 

IC0DE=1  264 

DEVELOP  U(N,M,K)  ARRAY  265 

DO  100  1=1,5  266 

IF  (MP.EO.NPOINT(I))  GO  TO  110  267 

100  CONTINUE  268 

PRINT  105, MP  269 

105  FORMAT  ( 'OINTEGRATION  CONTROL  =',I9,'  IS  NOT  AVAILABLE')  270 

STOP  271 

110  JF=KEY(I)  272 

JL=KEY(I+1)-1  273 

DO  115  1=1,5  274 

IF  (MPl.EO.NPOINT(I))  GO  TO  120  275 

115  CONTINUE  276 

PRINT  105,MP1  277 

STOP  278 

120  JF1=KEY(I)  279 

JL1=KEY(I+1)-1  280 

PD2=PIE/2.D0  281 
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K=1  282 

DO  130  J=JF1,JL1  283 

IF  (Y(J).NE.O.DO)  GO  TO  125  284 

THET1(K)=PD2  285 

WTTH(K)=WT(J)  286 

SINTH1{K)=1.D0  287 

C0STH1(K)-.0.D0  288 

GO  TO  130  289 

125  PDY=PD2*YiJ;  290 

THET1(K)=P02+PDY  291 

WTTH(K)=WT(J)  292 

SINTH1(K)=DSIN(THET1(K))  293 

C0STH1(K)=DC0S(THET1(K))  294 

K=K+1  295 

THET1(K)=PD2-»PDY  296 

WTTH(K)=WT(J)  297 

SINTH1(K)=DSIN(THET1(K))  298 

C0STH1(K)=DC0S(THET1(K))  299 

130  K=K+1  300 

EOSQ=EO*EO  301 

MAX:=NMIN+15  302 

IF  (MAX. GT. 100)  MAX=100  303 

CLEAR  STORAGE  FOR  INTEGRALS  304 

DO  135  NSBF=1,NMAX  305 

DO  135  M=l,2  306 

DO  135  NRT=1,KMAX  307 

DEN(NSBF,M,NRT)=0.D0  308 

135  XNUM(NSBF,M,NRT)«0.00  309 

SUM  OVER  REGIONS  310 

DO  210  NREG=1,N0RG  311 

IREG=NREG  312 

PRECALCULATE  THE  POWER  DENSITY  TIMES  PHI  INTEGRAL  313 

FIRST  CALCULATE  RADIUS  DEPENDENT  PART  OF  SOURCE  TERM  314 

NN=(NREG-1)*NMIN  315 

R1=--S8DP(NREG)  316 

R2=0.DC  317 

IF  (NREG.GT.l)  R2=SBDP(NREG-1 ) 318 

R11=R1+R2  319 

R13=R1-R2  320 

RAVG=R13/2.D0  321 

RCP=RHOP(NREG)*CP(NREG)  322 

FAC=EOSQ*. 5D0*S IGP (NREG)/4186 .03  323 

J=0  324 

DO  180  J3=JF,JL  325 

IF  (Y(J3).NE.0.D0)  GO  TO  140  326 

13=1  327 

ARG3(1)=R11  328 

GO  TO  145  329 

140  13=2  330 

R12=R13*Y(J3)  331 

ARG3(1)=R11+R12  332 

ARG3(2)=Rli-R12  > 333 

145  DO  180  L=1,I3  334 

R=ARG3(L)/2.D0  335 

J=J+1  336 

RR(J)=R  337 

R=R/100.D0  338 
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Q=R*FKP(NREG)  339 

CALL  BJYH{BJNP,BHNP,Q,NC,ST0PR,NMIN+2)  340 

NC=MIN0(NC-2,NMIN)  341 

NCK=0  342 

DO  150  N=1,NC  343 

FAC1=2*N+1  344 

NNN=NN+N  345 

W=BNP(NNN)  346 

X=BJNP{N+1)  347 

Y1=BETNP(NNN)  348 

Z=BHNP(N+1)  349 

NCK=NCK+1  350 

‘ T=FAC1*(W*X+Y1*Z)  351 

CALL  TERM(NCK,T,1)  352 

TR(N)=T  353 

T=ANP(NNN)*X+ALPNP(NNN)*Z  354 

CALL  TERM(NCK,T,0)  355 

TE(N)=r  356 

A=N+1  357 

B=N  358 

T=(W*fA*R.]NP(N)-B*BJNP(N+2))+Yl*(A*BHNP(N)-B^BHNP(N+2)))/FACl  359 

CALL  fERM(NU,T,l)  360 

TE‘(N)=T  361 

150  IF  (NCK.EQ.4)NCK=0  362 

THEN  CALCULATE  THETA  DEPENDENT  PART  OF  SOURCE  TERM  363 

DO  180  J2=1,MP1  364 

THETA=THET1(J2)  365 

SINTH=SINTH1(J2)  366 

C0STH=C0STH1(J2)  367 

CALL  PL  368 

ERAD=DCMPLX(O.DO,O.DO)  369 

ETHETA=DCMPLX(O.DO,O.DO)  370 

EPHI=DCMPLX{O.DO,O.DO)  371 

ITP=0  372 

TR=n  ‘^7'^ 

DO  165  N=1.NC  374 

FAC1=2*N+1  375 

IF  (IR.EQ.l)  60  TO  155  376 

T=P(N)*TR(N)  377 

ERAD=ERAD+T  •—  378 

IF  (CDASS(T).LT.CDABS(ERAD)*EPS)IR=1  379 

155  IF  (ITP.EQ.l)  GO  TO  160  380 

NP1=N+1  381 

RATI0=FAC1/(N*NP1)  382 

A=RATIO*P{N)/SINTH  383 

B=-RATI0*DP(N)  384 

C=A*TE(N)+B*TE1(N)  385 

ETHETA=ETHETA+C  386 

T=A*TE1(N)+B*TE(N)  387 

EPHI=EPHI+T  388 

IF  (CDABS(C).LT.rDABS(ETHETA)*EPS.AND.COABS(T).LT.CDABS(EPHI)*EPS)  389 
1ITP=1  390 

160  IF  (IR+ITP.E0.2)  GO  TO  175  391 

165  CONTINUE  392 

PRINT  170,NMIN,NC,THETA,R,ST0PR,EPS  393 

170  FORMAT  {15X,'NMIN  =‘,13,'  NC  =',13,'  THETA  =',F9.6,'  R =',2PF9.6,'  394 

1 STOPR  =',IPD9.2,'  IS  TOO  SMALL  FOR  ACCURACY  OF', 09. 2)  395 


175  ERAD=ERAD/0  396 

STORE  SOURCE  TERM  TIMES  PHI  INTEGRAL  397 

ERT=DREAL{ERAD*DCONJG(FRAD)+ETHETA*DCONJG(ETHETA) ) 398 

EPl=DREAL(EPHI*DCONJG(EPHI))  399 

S(J,J2,1  )=FAC*PIE*(ERT+EP1)  ■ 400 

180  S(  J,J2,2  )=FAC*Pn2*(ERT-EPl)  401 

CALCULATE  NUMERATOR  AND  DENOMINATOR  INTEGRALS  402 

DO  210  NSBF=1,NMAX  403 

N1=NSBF-1  404 

DO  210  M=1,NSBF  405 

IF  (M.NE.1.AND.M.NE.3)  GO  TO  210  406 

M2=M/2+l  407 

M1=M-1  408 

DO  185  J=1,MP1  409 

185  ALP0L(J)=ALP(N1,M1,C0STH1(J))*SINTH1(J)  410 

DO  195  J3  = 1,MP  411 

SUMJ2=  O.DO  412 

DO  190  J2  = 1,MP1  413 

190  SUMJ2  = SUMJ2  + WTTH(J2)*S(J3,J2,M2)*ALP0L{J2)  414 

195  SUM2(J3,M2)  = SUMJ2  415 

SUM2  IS  THE  INTEGRAL  OF  THE  SOURCE  TERM  TIMES  ALPOL  416 

ALPOL  IS  THE  PRODUCT  OF  THE  LEGENDRE  POLYNOMIAL  TIMES  THE  417 

SINE  OF  THETA  418 

J2  IS  AN  INDEX  FOR  THE  THETA  COORDINATE  ASSOCIATED  WITH  419 

GAUSSIAN  INTEGRATION  420 

DO  205  NRT=1,KMAX  421 

INTEGRATE  OVER  RADIUS  422 

F=FACT(IREG,NRT,NSBF)  423 

S1=XLAMDA(NRT,NSBF)*RCP-BP(IREG)  424 

SUM=O.DO  425 

'',UM3=0.D0  426 

DO  200  j3=1,MP  427 

R=RR(J3)  428 

R1=R  429 

CALL  SRBF(XJ,XY,DJ,DY)  430 

ZZ=AJ(NREG,NRT,NSBF)*XJ  +BY(NREG,NRT.NSBF)*XY  431 

IF  (DABS{XY).GT.1,D34)  PRINT  666,NSBF,M,NRT,AJ (NREG,NRT,NSBF),  432 

1 XJ,BY(NREG,NRT.NSBF),XY,ZZ  433 

666  FORMAT  (3I5,1P7D15.7)  434 

RSO=R*R  435 

INTEGRATE  OVER  THETA  436 

J=JF+(J3-i)/2  437 

WTJ=WT(J)*ZZ*RSQ  438 

SUM=SUM+WTJ*ZZ  439 

200  SUM3  = SUMS  + WTJ*SUM2(J3,M2)  440 

DEN(KSBF.M2,NRT)=DEN(NSBF,M2,NRT)+RCP*SUM*RAVG  441 

205  XNUM(NSBF,M2,NRT)=XNUM(NSBF,M2.NRT)+SUM3*RAVG  442 

210  CONTINUE  443 

CALCULATE  COEFFICIENTS  U(N,M,K)  444 

215  IF  (TCUT.GT.TIME)  TCUT=TIME  445 

IA=TCUT/TBPER  446 

IC=(TCUT-IA*TBPER)/TPER  447 

IB=MINO(NPUL,IC)  448 

XL=IA*TBPER+I8*TPER  449 

D=O.DO  450 

IF  {IC.LT,NPUL)  D=l.  451 

XU=DMIN1(TCUT,XL+D*TDUR)  452 


TA=TIME-T0UR-(IA-1)*TBPER-{NPUL-1)*TPER 

453 

TB=TIME-TDUR-IA*TBPER-(IB-1)*TPER 

454 

PRINT  220 

455 

220  FORMAT  ('OU  COEFFICIENTS') 

456 

DO  270  NSBF=1,NMAX 

457 

PRINT  225 

458 

225  FORMAT  ( ' ' ) 

459 

N1=NSBF-1 

460 

DO  270  M=1.NSBF 

461 

IF  (M.NE.1.AN0.M.NE.3)  GO  TO  270 

462 

M1=M-1 

463 

M2=M/2+l 

464 

NMM=N1-M1 

465 

NPM=N1+M1 

466 

F=1.D0 

467 

IF  (Ml.EO.O)  GO  TO  255 

468 

IF  (Nl.NE.Ml)  GO  TO  240 

469 

DO  235  I=2»NPM 

470 

235  F=F*I 

471 

GO  TO  250 

472 

240  II=2*M1 

473 

F1=NMM+1 

474 

DO  245  1=1,11 

475 

F=F*F1 

476 

245  F1=F1+1.D0 

477 

250  F=1.D0/F 

478 

255  F=(2.D0*N1+1.00)/(2.00*PIE)*F*PD2 

479 

DO  260  NRT=1,KMAX 

480 

IF  (M.NE.l)  GO  TO  260 

481 

XR=XLAMOA(NRT,NSBF) 

482 

D=0.D0 

483 

IF  (lA+IB.EQ.O)  GO  TO  258 

484 

X1=XR*TPER 

485 

X3=1.D0 

486 

IF  (X1.LE.40,D0)  X3=1.D0"DEXP(-X1) 

487 

IF  (lA.LE.O)  GO  TO  256 

488 

X4=XR*TA 

489 

IF  {X4.GT.87.D0)  GO  TO  256 

490 

X1=XR*NPUL*TPER 

491 

X5=1.D0 

492 

IF  {X1.LE.40.D0)  X5=1.D0-0EXP{-X1) 

493 

X1=XR*TBPER 

494 

X6=1.D0 

495 

X7=1.D0 

496 

IF  (Xl.GT.40.00)  60  TO  261 

497 

X6=1.00-DEXP(-X1) 

498 

X1=X1*IA 

499 

IF  (X1.LE.40.D0)  X7=l.nO-OEXP(-Xl) 

500 

261  0=D+DEXP(-X4)'X5*X7/{X3*X6) 

501 

256  IF  (IB.LE.O)  GO  TO  257 

502 

X4=XR*T8 

503 

IF  (X4.GT.87.00)  GO  TO  257 

504 

X1=XR*IB*TPER 

505 

X5=1.00 

506 

IF  {Xi.LE.40.no)  X5=1.D0-DEXP(-.X1) 

507 

0=D+DEXP{-X4)<»X5/X3 

508 

257  X1=XR»TDUR 

509 
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X4=1.D0  510 

IF  (XI. LE. 40. DO)  X4-1.D0-DEXP(-X1)  511 

D=D*X4  512 

258  IK  (XU.LE.XL)  GO  TO  259  513 

X1=XR*(TIMF.-XU)  514 

IF  (XI. GE. 87. DO)  GO  TO  259  515 

X3=XR*(XU-.XL)  516 

X4=1.D0  517 

IF  (X3.LE.40.D0)  X4=1.D0-DEXP(-X3)  518 

D=0+DEXP(-X1)*X4  519 

259  ETIME(NRT)=D/XR  520 

260  U (NSBF ,M2 ,NRT )=ETIME (NRT) ♦F*XNUH(NSBF ,M2 ,NRT)/DEN (NSBF ,M2 ,NRT)  521 

PRINT  265,N1,M1,(U(NSBF,M2,K),K=1,KMAX)  522 

265  FORMAT  (2I3,1P10D12.4/(8X,10D12.4))  523 

270  CONTINUE  524 

***  ABSORBED-POWER  DENSITY  AND  TEMPERATURE  RISE  AT  525 

***  GIVEN  POINTS  INTERIOR  TO  P-TH  REGION  526 

IF  (ISAR.NE.O)  ISAR=1  527 

PRINT  275,ZLAB(ISAR+1)  528 

275  FORMAT  ( 'O' ,29X, ' INTERNAL  POINT' .IIX, 'ABSORBED  POWER ',7X, 'TEMPERAT  529 
1URE',8X, 'APPROXIMATE'  530 

1 /IIX, 'POINT  REGION  RADIUS  THETA  PHI' .12X, 'DENSITY' ,14X,  531 
l'RISE',14X, 'ERROR'  532 

1/28X,'CM  DEG  DEG' ,12X,A7,13X. 'DEG  C'.13X,'PER  CENT'/)  533 
DO  345  II=1,N0CR  534 

READ  (5,30)  R,THETAD,PHID  535 

R R-COORDINATE  OF  PT  536 

ThFTAD  THETA  COORDINATE{DEGREES)  537 

PHID  PHI-COORDINATE (IN  EQUATORIAL  PLANE) (DEGREES)  538 

IF  (R.LE.O.DO)  GO  TO  290  539 

DO  285  Nf;EG=l,NORG  540 

IF  (R.LE.SBDP(NREG))  GO  TO  300  541 

285  CONTINUE  542 

290  NREG=1000000000  543 

PRINT  295,II,NREG,R,THETAD,PHID  544 

295  FORMAT  (I14,I8,3F10.3, ' **  THE  RADIUS  IS  OUTSIDE  THE  SPHERE  **')  545 

GO  TO  345  546 

NREG  = NUMBER  OF  THE  REGION  IN  WHICH  TEMP  IS  TO  BE  COMPUTED  547 
300  IREG=NREG  548 

R1=R  549 

R=R/1.D2  55C 

THETA=THETAD/RAD  551 

PHI=PHID/RAD  552 

CALL  BJYH(BJNP,BHNP,FKP(NREG)*R,NC,ST0PR,NMIN+2)  553 

NC=MIN0(NC-2,NMIN)  554 

SINTH=DSIN(THETA)  55£ 

COSTH=DCOS (THETA)  556 

CALL  PL  557 

CALL  EVEC(PD)  556 

PD=.5D0*SIGP(NREG)*PD  559 

KMAX1=KMAX  56C 

K1=KMAX-1  56] 

DO  315  KMAX=K1,KMAX1  562 

TRM=O.DO  563 

DO  315  NSRF=1,NMAX  56^ 

N1=NSBF-1  56E 

DO  315  M=1,NSBF  566 


IF  (M.NE.1.AND.M.NE.3)  GC  VO  315  567 

568 

M2=M/2+l  569 

ALPNM=ALP(N1 ,M1 ,COSTH)*OCOS(Hl*PHI ) 570 

IF  (ALPNM.EO.O.DO)  GO  TO  310  571 

SUM=O.DO  572 

DO  305  NRT=1,KHAX  573 

S1=XLAM0A(NRT,NSBF)*RH0P(IREG)*CP{IREG)-BP(IREG)  574 

F=FACT(IRE6,NRT,NS8F)  575 

CALL  SRBF(XJ,XY,OJ,DY)  576 

305  SUM=SUM+U(NSBF,M2,NRT)*{AJ(NRE6,NRT,NSBF)*XJ+BY{NREG,NRT,NSBF)*XY)  577 
310  TRM=TRM+SUM*ALPNM  578 

IF  (M.NE.3)  GO  TO  315  579 

IF  (KMAX.EQ.Kl.AND.NSBF.EO.NMAX)  SRK1*TRM  580 

IF  (KMAX.EQ.KMAXl.AND.NSBF.EQ.NMAX-l)  SBFMl^TRM  581 

315  CONTINUE  582 

KMAX=KMAX1  583 

PCER=(TRM-SRM1)/TRM  584 

PCE8F=(TRM-SBFM1)/TRM  585 

IF  (ISAR.E0.O)PD=PD/RHOP(NREG)  586 

♦**  PRINT  PARTICULARS  OF  INTERIOR  POINT  OF  REGION  P 587 

PRINT  340,II.NREG,R1.THETAD,PHI0.PD,TRM.PCEBF.PCER  588 

340  FORMAT  (I14,I8.F10.3,2F8.2,F19.8.1PD20.4,2P2F14.7)  589 

345  CONTINUE  590 

IF  {IPL1.E0.0.AND.IPL2.EQ.0)  GO  TO  10  591 

IF  (IPLSW.EQ.l)  GO  TO  350  592 

IPLSW=1  593 

CALL  PL0TS(0,0,8)  594 

CALL  PL0T(0.,-ll.,-3)  595 

CALL  PLOT (0., 2., *3)  596 

NTR=MAX0(NTR,1)  597 

NTR=MIH0{NTR,5)  598 

350  IF  (IPLl.EO.O)  GO  TO  405  599 

***  PLOT  POWER  DENSITIES  ALONG  Z,  X AND/OR  Y AXIS  600 

NPTS=300  601 

NPT02=NPTS/2  602 

NP2=NPT02+1  603 

OX=SBDP(N0RG)/NPT02  604 

DO  400  KJI=1,3  605 

IF  (KJ!.E0.1.ANO.(IPL1.E0.2.0R.IPL1.E0.3.0R.IPL1.E0.6))  GO  TO  400  606 

IF  (KJI.EQ.2.AND.(IPL1.EQ.1.0R.IPL1.EQ.3.0R.IPL1.EQ.5))  GO  TO  400  607 

IF  (KJI.E0.3.AND.(IPL1.EQ.1.0R.IPL1.E0.2.0R.IPL1.E0.4))  60  TO  400  608 

PRINT  355  609 

355  FORMAT  (’O')  610 

TRMAX=0.  611 

C0STH=0.00  612 

IF  (KJI.EO.l)  C0STH=1.D0  613 

IRE6=N0RG  614 

R1=SBDP(N0RG)  615 

***  CALCULATE  POWER  DENSITIES  ALONG  SPHERE  DIAMETER  616 

DO  370  I=1,NP2  617 

RC=RHOP(IREG)*CP{IREG)  618 

BP1=BP{IREG)  619 

TRM=O.DO  620 

TRM1=0.D0  621 

DO  365  NSBF=1,NMAX  622 

N1=NSBF-1  623 
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C0SMP=1.D0  624 

DO  365  M=1,NSBF  625 

IF  (M.NE.1.AND.M.NE„3)  GO  TO  365  626 

IF  (KJI.E0.3.AND.M.EQ.3)  COSMP=-l.DO  627 

M1=M-1  628 

M2=M/2+l  629 

SUM=O.DO  630 

DO  360  NRT»1,KMAX  631 

S1=XLAMDA(NRT.NSBF)*RC-BP1  632 

F=FACT(IREG,NRT,NSBF)  633 

CALL  SRBF(XJ,XY.DJ,DY)  634 

360  SUM=SUM+U(NSBF,M2,NRT)*(AJ(IREG.NRT.NSBF)*XJ+BY(IREG,NRT,NSBF)*XY)  635 
TRM=TRM+SUM*ALP{N1 ,M1 ,C0STH)*C0SMP  636 

TRM1=TRM1+SUM*ALP(N1 ,M1 ,-COSTH)*COSMP  637 

365  CONTINUE  638 

R3(I)=R1  639 

TR3(I)=TRM  640 

R3(NPTS-I+3)=-Rl  641 

TR3(NPTS-I+3)=TRM1  642 

TRMAX=DMAX1(TRM,TRM1,TRMAX)  643 

R1=R1-DX  644 

IF  {IRE6.GT.1.AND.R1.LT.SBDP(IREG-1))IREG=IREG-1  645 

370  IF  (R1.LT..0001)R1=.0001  646 

***  DETERMINE  PLOT  SCALE  FOR  POWER  DENSITIES  647 

PD3=,0001  648 

DO  375  1=1,10  649 

PD3=5.*PD3  650 

IF  (TRMAX.LT.PD3)  60  TO  380  651 

PD3=PD3*2.  652 

IF  (TRMAX.LT.PD3)  GO  TO  380  653 

375  CONTINUE  654 

380  TRMAX=PD3  655 

***  PLOT  POWER  DENSITY  ALONG  DIAMETER  ON  Z,  X OR  Y AXIS  656 

BLAB(1)=AX(KJI)  657 

DO  390  1=1, NTR  658 

ANG=2*(I-1)*PIE/NTR  659 

AX1=.01*COS(ANG)  660 

AY=.01*SIN{ANG)  661 

IF  (NTR.EO.l)  AX1=0.  662 

CALL  PlOT(AXl,AY,-3)  663 

390  CALL  PLTCV1(R3,TR3,5.,6.,BLAB,DLAB,22,26,NPTS+2,0,1,1,-R3(1),  664 

1R3(1),0.,TRMAX,0,0,.14,R3(1)/3.,TRMAX/5.,1)  665 

CALL  PL0T(7.,n.,-3)  666 

400  CONTINUE  667 

405  IF  (IPL2.EQ.0)  GO  TO  10  668 

***  PLOT  POWER  DENSITY  CONTOURS  IN  E PLANE,  H PLANE  AND/OR  X-Y  PLANE  669 

NPTS=100  670 

NPTD2=NPTS/2  671 

NPTP2=NPTS+2  672 

X1=SBDP(N0RG)  673 

XF=10./{2.*X1)  674 

DX=X1/NPTD2  675 

X3=X1  676 

DO  410  I=1,NPTD2  677 

X2(I)=X3  678 

X2(NPTS+3-I)=-X3  679 

410  X3=X3-0X  680 


X2{NPTD2+1)=.0001  681 

X2(NPTD2+2)=-.0001  682 

***  CALCULATE  POWER  DENSITIES  AT  POINTS  IN  PLANE  683 

N12=NPTP2/2  684 

DO  465  KJI=1,3  685 

IF  (KJI.E0.1.AND.(IPL2.E0.2.0R.IPL2.E0.3.0R.IPL2.E0.6))  GO  TO  465  686 

IF  (KJI.E0.2.AND.{IPL2.EQ.1.0R.IPL2.E0.3.0R.IPL2.E0.5))  60  TO  465  687 

IF  (KJI.E0.3.AND.(IPL2.EQ.1.0R.IPL2.E0.2.0R.IPL2.EQ.4))  60  TO  465  688 

Y3=0.  689 

X3=0.  690 

Z3=0.  691 

DO  455  1=1, N12  692 

Il=NPTS+3-I  693 

DO  455  J=l,N12  694 

Jl=NPTS+3-J  695 

IF  (KJI.GT.l)  GO  TO  415  696 

X3=X2(I)  697 

Z3=X2(J)  698 

GO  TO  425  699 

415  IF  (KJI.GT.2)  GO  TO  420  700 

Y3=X2(I)  701 

Z3=X2(J)  702 

GO  TO  425  703 

420  X3=X2(I)  704 

Y3=X2(J)  705 

425  R1=DSQRT(X3*X3+Y3*Y3+Z3*Z3)  706 

IF  IRI.LE.XI)  60  TO  430  707 

DAR(I,J)=-1.  708 

DAR(I,J1)=-1.  709 

GO  TO  453  710 

430  DO  435  IREG=1,N0RG  711 

IF  (Rl.LE.SBDP(IREG))  GO  TO  440  712 

435  CONTINUE  713 

***  CALCULATE  TEMPERATURE  RISE  AT  POINTS  IN  PLANE  714 

440  C0STH=Z3/R1  715 

PHI=DATAN2(Y3,X3)  716 

RC=RHOP(IREG)*CP(IREG)  717 

BP1=BP(IRE6)  718 

TRM=O.DO  719 

TRM1=0.D0  720 

DO  450  NSBF=1,NMAX  721 

N1=NSBF-1  722 

DO  450  M=1,NSBF  723 

IF  (K.NE.1.AND.M.NE.3)  GO  TO  450  724 

M1=M-1  725 

M2=M/2+l  726 

SUM=O.DO  727 

DO  445  NRT=1,KMAX  728 

S1=XLAMDA(NRT,NSBF)*RC-BP1  729 

F=FACT(IREG,NRT,NSBF)  730 

CALL  SRBF(XJ,XY,DJ,DY)  731 

445  SUM=SUM+U(NSBF,M2,NRT)*(AJ(IRE6,NRT,NSBF)*XJ+BY(IRE6,NRT,NSBF)*XY)  732 
TRM=TRM+SUM*ALP{N1 ,M1 ,COSTH)*OCOS(M1*PHI ) 733 

TRM1=TRM1+SUH*ALP(N1,M1,-C0STH)*0C0S(M1*PHI)  734 

450  CONTINUE  735 

DAR(I,J)=TRM  736 

DAR(I,J1)=TRM1  737 


126 


453  DAR(I1,J)=DAR(I,J)  738 

455  DAR(I1,J1)=DAR(I,J1)  739 

***  PLOT  CONTOURS  740 

no  460  1=1, NTR  741 

ANG=2*(I-1)*PIE/NTR  742 

aX1=.01*C0S(ANG)  743 

AY=.01*SIN(ANG)  744 

IF  (NTR.EQ.l)  AX1=0.  745 

CALL  PL0T(AXl.AY,-3)  >'% 

CALL  SYMB0L{-.5,6...21,CLAB(1,KJI),0..9)  747 

460  CALL  CNTRP1(X2.NPTP2,X2,NPTP2,DAR.10.0.IFL)  748 

CALL  PL0T(10.,0.,-3)  749 

465  CONTINUE  750 

GO  TO  10  751 

495  IF  (IPLSW.NE.O)  CALL  PL0T(0.,0..999)  752 

500  STOP  753 

END  754 

755 

756 

SUBROUTINE  COEF  757 

IMPLICIT  REAL *8  (A-H,0-Z)  758 

GENERATES  EXPANSION  COEFFICIENTS  759 

C0MPLEXM6  FKP,ANP,BNP,ALPNP,BETNP,BJNP.BHNP,BJHP1(500),BJHP2(500)  760 
1 .SJNPl (100) ,SHNP1 (100) ,DELNP.SNT11,SNT12,SNT21 .SNT22.TNT11 .TNT12,T  761 
1NT21,TNT22,ETAP1,ETAP2,ZEP1,ZEP2.SNP11.SNP12,SNP21,SNP22,TNP11,TNP  762 
112,TNP21,TNP22,0ELl,nEL2,RATI0,Z  763 

COMMON  FKP(7),ANP(300),BNP(300),ALPNP(300),BETNP(300),BJNP(100),BH  764 
1NP(100),BDP(6),P(51),DP(50).R.THETA.C0STH,PHI.SINTH.ST0PR,E0  765 

COMMON  /A/NORG,NREG,NRT,NSBF,NMIN,NC,ICODE  766 

DIMENSION  NTER(6)  767 

COMPUTE  EXPANSION  COEFFICIENTS  AN1,BN1,ANN,BNN,ALPN1,BETN1,  768 
ALPNN.BETNN  769 

N^=l  770 

NMIN=100  771 

DO  10  NR=1,N0RG  772 

CALL  BJYH(SJNP1,SHNP1,FKP(NR)*BDP(NR),N,ST0PR,NMIN)  773 

CALL  BJYH(BJNP,BHNP,FKP(NR+1)*B0P(NR),NN,ST0PR.NMIN)  774 

NMIN=MINO(N,NN,NMIN)  775 

N2=N1+NMIN  776 

DO  5 I=1,NMIN  777 

BJHP1(N1)=SJNP1(I)  778 

BJHPl(N2)=SHNPlfI)  779 

BJHP2(Nl)=BJNP(i)  780 

BJHP2(N2)=BHNP(n  781 

N1=N1+1  782 

5 N2=N2+1  783 

N1=N1+NMTN  784 

10  NTER(NR)=NMIN  785 

NMIi\’=NMIN-2  786 

IF  (NMIN.LE.50.ANn.N2.LE.301)  GO  TO  20  787 

PRINT  15,N2.NMIN  788 

15  FORMAT  ('OCOEF  ERROR:  N2  =’,13,'  NMIN  «’,I3,'  DIMENSIONS  ARE  TOO  789 

1 SMALL')  790 

STOP  791 

20  DO  25  1=1, NMIN  792 

ALPNP(I)=DCMPLX(O.DO,O.DO)  793 

25  BETNP(I)=DCMPLX(0.D0,0.D0)  794 
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NSUM=NORG*NMIN 

no  35  I=1,NMIN 

JJ=0 

KK=0 

XI  = I 

XI1=I+1 

XI2=2*I+1 

SNTll=DCMPLX(l.DO,O.DO) 

SNT12=DCMPLX(0.D0,0.D0) 

SNT21=SNT12 

SNT22=SNT11 

TNT11=SNT11 

TNT12=SNT12 

TNT21=SNT12 

TNT22=SNT11 

DO  30  J=1,N0RG 

KK=KK+NTER(J) 

KKI=KK+I 

JJI=JJ+I 

ETAP 1= ( X I 1 *B JHP 1 ( J J I ) -X I *B JHP 1 ( J J 1+2 ) ) /X I 2 

ETAP2={XIl*BJHP2(JJI)-XI*BJHP2(JJI+2))/XI2 

ZEPl=(XIl*BJHPl(KKI)-XI*BJHPl(KKI+2))/XI2 

ZEP2=(XIl*BJHP2(KKI)-XI*BJHP2(KKI+2))/XI2 

DELNP=BJHP1(JJI+1)*ZEP1-BJHP1(KKI+1)*ETAP1 

RATI0=FKP(J+1)/FKP(J) 

Z=RATI0+ETAP2 

SNPn=(ZEPl*BJHP2(JJI+l)-Z*BJHPl(KKI+l))/DELNP 
SNP21= (Z*BJHP1 (JJI+1 )-ETAPl*8JHP2(JJI+l ) )/DELNP 
Z*RATI0*ZEP2 

SNP12={ZEP1*BJHP2(KKI+1 )-Z*8JHPl (KKI+1 ) )/OELNP 
SNP22= (Z*BJHP1 (JJI+1 )-ETAPl*BJHP2(KKI+l ) )/DELNP 
Z=SNT11 

SNT11=SNT11*SNP11+SNT12*SNP21 

SNT12=Z*SNP12+SNT12*SNP22 

Z=SNT21 

SNT21=SNT21*SNP11+SNT22*SNP21 

SNT22=Z*SNP12+SNT22*SNP22 

Z=RATI0*ZEP1 

TNP11=(Z*BJHP2{JJI+1)-BJHP1(KKI+1)*ETAP2)/DELNP 
TNPl 2= (Z*BJHP2(KKI+1 )-BJHPl (KKI+1 )*ZEP2)/DELNP 
Z=RATI0*ETAP1 

TNP21=(BJHP1 (JJI+1 )*ETAP2-Z*BJHP2(JJI+1))/DELNP 
TNP22=(BJHP1 (JJI+1 )*ZEP2-Z*BJHP2(KKI+1))/DELNP 
Z=TNT11 

TNT11=TNT11*TNP11+TNT12*TNP21 

TNT12=Z*TNP12+TNT12*TNP22 

Z=TNT21 

TNT21=TNT21*TNP11+TNT22*TNP21 

TNT22=Z*TNP12+TNT22*TNP22 

JJ=JJ+2*NTER(J) 

30  KK=KK+NTER(J) 

ANP ( I ) =SNT1 1- (SNT12*SNT21 ) /SNT22 
BNP ( I )=TNT1 1- (TNT12*TNT21 )/TNT22 
LL=NSUM+I 

ANP(LL)=DCMPLX(1.D0.0.D0) 

BNP(LL)=DCMPLX(1.D0,0.D0) 

ALPNP(LL)=-SNT21/SNT22 


35  BETNP(LL)=-TNT21/TNT22  852 

IF  (NORG.EQ.l)  RETURN  853 

COMPUTE  EXPANSION  COEFFICIENTS  AN2,. . . .AN(N-l) ;BN2 8N(N-1  854 

);ALPN2 ALPN(N-1);BETN2,...,BETN(N-1)  855 

JJ=0  856 

KK=0  857 

MM1=0  858 

MM2=NMIN  859 

NRGM1=N0RG-1  860 

DO  45  J=1,NRGM1  861 

KK=KK+NTER(J)  862 

DO  40  I=1,NMIN  863 

KKI=KK+I  864 

JJI=JJ+I  865 

XI=I  866 

XI1=I+1  867 

XI2=2*I+1  868 

ETAPl={XIl*BJHPl(JJI)-XI*BJHPl(JJI+2))/XI2  869 

ETAP2=(XIl*BJHP2(JJI)-XI*BJHP2{JJI+2))/XI2  870 

ZEPl=(XIl*BJHPl(KKI)-XI*BJHPl(KKI+2))/XI2  871 

ZEP2=(XIl*BJHP2(KKI)-XI*BJHP2(KKI+2))/XI2  872 

DELNP=BJHP1(JJI+1)*ZEP1-BJHP1(KKI+1)*ETAP1  873 

RATI0=FKP(J+1)/FKP(J)  874 

Z=RATI0*ETAP2  875 

SNP11=(ZEP1*BJHP2(JJI+1)-Z*BJHP1(KKI+1))/DELNP  876 

SNP21=(Z*BJHP1 (JJI+1 )-ETAPl*BJHP2(JdI+l ) )/DELNP  877 

Z=RATI0*ZEP2  878 

SNP12=(ZEPl*BJHP2(KKI+l)-Z*BJHPl(KKI+l))/nELNP  879 

SNP22=(Z*BJHPl(JJI+l)-ETAPl*BdHP2(KKI+l))/0ELNP  880 

0EL1=SNP11*SNP22-SNP12*SNP21  881 

Z=RATI0*ZEP1  882 

TNP11=(Z*BJHP2(JJI+1)-BJHP1(KKI+1)*ETAP2)/DELNP  883 

TNP12={Z*BJHP2(KKI+1)-BJHP1(KKI+1)*ZEP2)/DELNP  884 

Z=RATI0*ETAP1  885 

TNP21=(BJHP1 (JJI+1 )*ETAP2-Z*BJHP2(JJI+1 ) )/DELNP  886 

TNP22=  (BJHPl  (JJI+1  )*ZEP2-Z*BJHP2(K!CI+1 ) )/DELNP  887 

DEL2=TNP11*TNP22-TNP12*TNP21  888 

NN1=MM1+I  889 

NN2=MM2+I  890 

ANP(NN2)=(ANP(NN1)*SNP22-ALPNP(NN1)*SNP12)/DEL1  891 

BNP (NN2 )= (BNP (NNl )*TNP22-BETNP (NNl )*TNP12 )/DEL2  892 

ALPNP (NN2 ) = (-ANP (NNl ) *SNP21+AlPNP (NNl ) *SNP 1 1 ) /DELI  893 

40  BETNP(NN2)=(-BNP(NNl)nNP21+BETNP(NNl)nNPll)/DEL2  894 

JJ=JJ+2*NTER(J)  895 

KK=KK+NTER(J)  ^ 896 

MM1=MM1+NMIN  897 

45  MM2=MM2+NMIN  898 

RETURN  899 

END  900 

901 

902 

SUBROUTINE  EVEC(PD)  903 

IMPLICIT  REAL *8  (A-H,0-Z)  904 

COMPUTES  THE  RADIAL, COLATITUDE.ANO  AZIMUTHAL  905 

COMPONENTS  OF  ELECTRIC  FIELD  VECTOR  E FOR  906 

REGION  P AND  SCALAR  PRODUCT  E.E*  907 

C0MPLEX*16  FKP,ANP,BNP,ALPNP,8ETNP,RJNP,BHNP,ERA0.ETHETA,EPHI ,T.T1  908 


1 C W X Y Z 909 

C0Mm6n’fKP(7),ANP(300),BNP(300),ALPNP(300),BETNP(300),BJNP(100),BH  910 
1NP(100),BDP(6),P(51),DP(50),R.THETA,C0STH,PHI,SINTH,ST0PR,E0  911 

COMMON  /A/NORG,NREG,NRT,NSBF,NMIN,NC,ICODE  912 

DATA  EPS/1. D-8/  913 

ERAD=DCMPLX(O.DO,0.00)  914 

ETHETA=DCMPLX(0.D0,0.DC)  915 

EPHI=DCMPLX(O.DO,O.DO)  916 

NCK=0  917 

NN=(NREG-1)*NMIN  918 

IR=0  919 

IF  (THETA. EQ.O. DO)  IR=1  920 

ITP=0  921 

DO  25  N=1,NC  922 

FAC1^2*N+1  923 

NNN=NN+N  924 

W=BNP(NNN)  925 

X=BJNP(N+1)  926 

Y=BETNP(NNN)  927 

Z=BHNP(N+1)  928 

NCK=NCK+1  929 

IF  (IR.EQ.l)  GO  TO  5 930 

T=FAC1*P(N)*(W*X+Y*Z)  931 

CALL  TERM(NCK,T,1)  932 

ERAO=ERAD+T  933 

IF  {CDABS(T).LT.CDABS(ERAn)*EPS)IR*l  934 

5 IF  (ITP.EQ.l)  GO  TO  20  935 

T=ANP(NNN)*X+ALPNP(NNN)*Z  936 

CALL  TERM{NCK,T,0)  937 

NP1=N+1  938 

RATI0=FAC1/(N*NP1)  939 

A=NP1  940 

B»N  941 

Tl= (W* (A*8 JNP (N )-B*BJNP (N+2 ) )+Y* (A*BHMP (N )-B*BHNP (N+2 ) ) ) /FACl  942 

CALL  TrRM(NCK,Tl,l)  943 

IF  (SINTH.GT.l.D-6)  GO  TO  10  944 

A-FAC1/2.D0  945 

IF  (THETA.GE.3.141591DQ)A=A*(-1.D0)*'»KP1  946 

GO  TO  15  947 

10  A=RATIO*P(N)/SINTH  948 

15  l-.RATIO*DP(N)  949 

C-A*T+B*T1  950 

ErHETA=ETHETA+C  951 

T=A*T1+B*1  952 

EPHI=EPHI+T  953 

IF  (CDABS{C).LT.CDA8S(ETHETA)*EPS.AND..CDABS(T).LT.CDABS(EPHI)*EPS)  954 
1ITP=1  955 

20  IF  {IR+ITP.E0.2)  GO  TO  35  956 

25  IF  (NCK.E0.4)NCK=0  957 

PRINT  30,NMIN,NC,ST0PR,EPS  958 

30  FORMAT  {15X,'NMIN  =',13,'  NC  =',I3,'  STOPR  =',1PD14.4,'  IS  TOO  SM  959 
lALL  FOR  ACCURACY  0F',D14.4)  960 

35  ECOSPH=EO*OCOS(PHI)  961 

ERAD=-ECOSPH/ (FKP (NREG) ♦R )*ERAD  962 

£THETA=ECOSPH*ETHETA  963 

EPHI=EO*OSIN{PHI)*EPHI  964 

PD=DREAL (ERAD*DCONJG(ERAD ) )+DREAL(ETHETA*DC0NJ6{ETHETA) )+OREAL (EPH  965 


1I*DC0NJG{EPHI))  966 

RETURN  967 

END  968 

969 

970 

SUBROUTINE  TERM(NCK,T,KEY)  971 

IMPLICIT  REAL*8  (A-H,0-Z)  972 

COMPUTES  I**NCK*(N-TH  TERM  IN  SERIES)  973 

C0MPLEX*16  T 974 

IF  (KEY.EQ.l)  GO  TO  5 975 

GO  TO  {10,15,20,25),NCK  976 

5 GO  TO  (15,20,25,10),NCK  977 

10  T=OCMPLX(-DIMAG(T),DREAL(T))  978 

GO  TO  25  979 

15  T=-T  980 

GO  TO  25  981 

20  T=DCMPLX(DIMAG(T),-DREAL(T))  982 

25  RETURN  983 

END  984 

985 

986 

SUBROUTINE  PL  987 

IMPLICIT  REAL*8  {A-H,0-Z)  988 

ASSOCIATED  LEGENDRE  FUNCTIONS  OF  THE  FIRST  989 

KIND,  DEGREE  K AND  ORDER  1 AND  THEIR  FIRST  990 

DERIVATIVES  991 

C0MPLEX*16  FKP,ANP,BNP,ALPNP,BETNP,BJNP,BHNP  992 

COMMON  FKP(7),ANP(300),BNP(300).ALPNP(300).BETNP(300).BJNP{100),BH  993 
lNP(100),BDP(6).P(!jl),DP(50),R,THETA,C0STH,PHI,SINTH,ST0PR,E0  994 

COMMON  /A/NORG,NREG,NRT.NSBF,NMIN,NC,ICOOE  995 

P(1)=SINTH  996 

P(2)=3.D0*SINTH*C0STH  997 

DP(1)=C0STH  998 

DO  10  M=2,NC  999 

A=M  1000 

MP1=M+1  1001 

P(MP1)=(2.D0*A+1.D0)/A*C0STH*P(M)-(A+1.D0)/A*P(M-1)  1002 

IF  (SINTH.GT.l.D-6)  GO  TO  5 1003 

DP(M)=M*MPl/2  1004 

IF  (THETA.GE.3.141591D0)DP{M)^{-1.D0)**M*DP(M)  1005 

GO  TO  10  1006 

5 DP(M)=(A*C0STH*P(M)-{A+1.D0)*P(M-1))/SINTH  1007 

10  CONTINUE  1008 

RETURN  1009 

END  1010 

1011 

1012 

FUNCTION  ALP(N,M,X)  1013 

IMPLICIT  REAL*8  (A-H,0-Z)  1014 

ASSOCIATED  LEGENDRE  FUNCTIONS  OF  THE  FIRST  KINO,  1015 

DEGREE  N AND  ORDER  M.  N AND  M GTE  0,  N GTE  M 1016 

FM=M  1017 

IF  (M.GT.O)  GO  TO  5 1018 

P1=1.D0  1019 

GO  TO  25  1020 

5 IF  (M.GT.l)  GO  TO  10  1021 

SUM=2.D0  1022 


m 

>V« 

^ 

S•!^ 

i 


GO  TO  20  1023 

10  J=2*M  1024 

SUM=J  1025 

IS=M-1  1026 

DO  15  1=1, IS  1027 

15  SUM=SUM*(J-I)  1028 

20  P1=SUM*(('1.D0-X*X)**(FM/2.D0))/(2.D0**M)  1029 

25  IF  (N.NE.M)  GO  TO  30  ’030 

ALP=P1  1031 

GO  TO  40  1032 

30  ALP=(2.D0*FM+1.DC)*X*P1  1033 

IF  (N.EQ.M+1)  GO  TO  40  1034 

IS=N-M  1035 

DO  35  1=2, IS  1036 

P2=ALP  1037 

C1=2*(M+I)-1  1038 

ALP={C1*X*P2-(C1-I)*P1)/I  1039 

35  P1=P2  1040 

40  RETURN  1041 

END  1042 

1043 

1044 

SUBROUTINE  RFNDR(RSTART,STEP1,E,NRTS,M1,NITR)  1045 

EXPLICIT  REAL*8  (A-H,0-Z)  1046 

ROOT  FINDER  1047 

COMMON  /A/NORG,NREG,NRT,NSBF,NMIN,NC,ICODE  1048 

COMMON  /B/FACT(6,25,18),AJ(6,25,18),[r(6,25,18),XLAMDA(25,18),SBDP  1049 
1(6),RH0P(6),CP(6),BP{6),TCP(6),H  1050 

EXTERNAL  FNCAL  1051 

STEP=STEP1  1052 

M=Ml-3  1053 

1=1  1054 

SL=RSTART  1055 

5 X=SL  1056 

NRT=I  1057 

W=FNCAL(X)  1058 

10  IF  (W)  15,55,25  1059 


15  DO  20  J=1,M 

1060 

X=X+STEP 

1061 

tv-' 

w"-- 

V=FNCAL(X) 

1062 

IF  (V)  20,55,50 

1063 

20  W=V 

1064 

GO  TO  35 

1065 

**• 

25  DO  30  J=1,M 

1066 

X=X+STEP 

1067 

V=FNCAL(X) 

1068 

IF  (V)  50,55,30 

1069 

30  W=V 

1070 

' c.*-* 

35  IF(M.GT.IOOO)  GO  TO  40 

1071 

M=M+1 

1072 

STEP=STEP*1.D1 

1073 

S.--' 

GO  TO  10 

1074 

''‘ti 

40  PRINT  45,SL,X 

1075 

45  FORMAT  ('ORFNDR  ERROR:  NO  ROOTS  FROM' ,1PE14.4, ' T0',E14.4) 

1076 

STOP 

1077 

50  SL=X-STEP 

1078 

SR=X 

1079 

132 


CALL  DRTMI(X,F,FNCAL,SL,SR,W,V.E.NITR) 

55  XLAMDA(I,NSBF)=X 
SL=X+STEP 

IF(I+NSRF.EQ.2)  STEP  =DMAX1(X/10.00,STEP) 
IF  (I.GT.1)STEP=(X-XLAMDA(I-1,NSBF))/10.D0 
1=1+1 

IF  (I.LE.NRTS)  GO  TO  5 

RETURN 

END 


SUBROUTINE  DRTMI (X.F.FCT.XLI ,XRI.FLI,FRI.EPS.IEND) 
IMPLICIT  REAL *8  (A-H.O-Z) 

BISECTION  METHOD 

XL=XLI 

XR=XRI 

FL=FLI 

FR=FRI 

1=0 

T0LF=100.D0*EPS 
5 1=1+1 

DO  30  K=1,IEND 
X=.5D0*(XL+XR) 

F=FCT(X) 

IF  (F.EQ.O.DO)  GO  TO  45 

IF  (DSIGN{1.D0,F)+DSIGN(1.D0,FR).NE.0.D0)  GO  TO  JO 
TOL=XL 
XL=XR 
XR=TOL 
TOL=FL 
FL=FR 
FR=TOL 
10  TOL=F-FL 
A=F*TOL 
A=A+A 

IF  (A.GE.FR*(FR-FL))  GO  TO  25 
IF  (I.GT.IEND)  GO  TO  25 
A=FR-F 

DX= ( X-XL )*FL* ( 1 .DO+F*(A-TOL ) /(A*(FR-FL ) ) ) /TOL 

XM=X 

FM=F 

X=XL-DX 

F=FCT(X) 

IF  (F.EQ.O.DO)  GO  TO  45 

TOL=EPS 

A=DABS(X) 

IF  (A.GT.1.D0)T0L=T0L*A 
IF  (DABS(DX).GT.TOL)  GO  TO  15 
IF  (DABS(F).LE.TOLF)  GO  TO  45 
15  IF  {DSIGN(1.D0,F)+DSI6N(1.00,FL).NE.0.D0)  GO  TO  20 
XR=X 
FR=F 
GO  TO  5 
20  XL=X 
FL=F 
XR=XM 
FR=FM 
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1084 

1085 

1086 

1087 
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1092 

1093 

1094 

1095 

1096 

1097 

1098 

1099 

1100 
1101 
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1103 

1104 

1105 

1106 

1107 

1108 

1109 
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1112 

1113 

1114 

1115 

1116 

1117 

1118 

1119 
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1121 
1122 

1123 

1124 

1125 

1126 

1127 

1128 

1129 

1130 

1131 

1132 

1133 

1134 

1135 
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GO  TO  5 1137 

25  XR=X  1138 

FR=F  1139 

TOL=EPS  1140 

A=DABS(XR)  1141 

IF  (A.GT.1.D0)T0L=T0L*A  1142 

IF  (DABS(XL-XR).GT.TOL)  GO  TO  30  1143 

IF  (DA8S(FR-FL).LE,T0LF)  GO  TO  40  1144 

30  CONTINUE  1145 

PRINT  35,XL,XR  1146 

35  FORMAT  ('ODRTMI  ERROR.'  ROOT  BETWEENMPD15.7, * AND',D15.7,'  MAY  BE  1147 
1 INACCURATE')  1148 

40  IF  {DABS(FR).LE.DABS(FL))  GO  TO  45  1149 

X=XL  1150 

F=FL  1151 

45  RETURN  1152 

END  1153 

1154 

1155 

FUNCTION  FNCAL(EIGV)  1156 

FUNCTION  EVALUATOR  USED  IN  THE  DETERMINATION  1157 

OF  THE  EIGENVALUES  LAMBDANK  1158 

IMPLICIT  REAL*8  (A-H,0-Z)  1159 

COMMON  /B/FACT(6,25,18),AJ(6,25,18),BY(6,25,18),XLAMDA(25,18),SBDP  1160 
1{6),RH0P(6).CP(6),BP(6),TCP(6),H  1161 

COMMON  /C/AJ1,S,F,R,I  1162 

COMMON  /A/NORG,NREG,NRT,NSBF,NMIN,NC,ICODE  1163 

BY(1,NRT,NSBF)  = O.DO  1164 

DO  35  I = 1,N0RG  1165 

S =(EIGV*RHOP(I)*CP(I)-BP(I))/TCP(I)  1166 

F=DSQRT(DABS(S))  1167 

FACT(I,NRT,NSBF)=F  1168 

IF  (I.NE.l)  60  TO  27  1169 

IF  (F.NE.O.DO)  GO  TO  5 1170 

AJ(1,NRT,NSBF)=AJ1  1171 

GO  TO  30  . 1172 

5 AJ(1,NRT,NSBF)-AJ1/F**(NSBF-1)  1173 

IF  (S.LT.O.DO)  AJ(l,NRT,NSBF)=AJ(l,NRT,NSBF)/((-l)**((NSBF>l)/2))  1174 
GO  TO  30  1175 

27  R=SBDP(I-1)  1176 

CALL  SRBF(AM,BM,ATM,BTM)  1177 

DELTA  = AM*BTM-ATM*BM  1178 

T1  = AJ(I-1,NRT,NSBF)*A  + BY(I-1.NRT,NSBF)*BE  1179 

T2  = AJ(I-1,NRT,NSBF)*AT  + BY(I-1,NRT,NSBF)*BT  1180 

AJ(I,NRT,NSBF)  = (T1*BTM-T2*BM) /DELTA  1181 

BY(I,NRT,NSBF)  = (T2*AM-T1*ATM) /DELTA  1182 

30  R=SBDP(I)  1183 

35  CALL  SRBF(A,BE,AT,BT)  1184 

FNCAL =A J ( NORG , NRT . N SBF ) *AT+B Y (N0R6 , NRT , NSBF ) *BT  1185 

1 +H*(AJ(N0RG,NRT,NSBF)*A+BY{N0R6,NRT.NSBF)*BE)  1186 

RETURN  1187 

END  1188 

1189 

1190 

SUBROUTINE  BJYH(BJNP,BI;NP,Z,N,ST0PR,NBF)  1191 

IMPLICIT  C0MPLEX*16{A-H,0-Z)  1192 

DIMENSION  BJNP(62),BHNP(62)  1193 
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REAL *8  STOPR,X»XNPH,DREAL,niMAG  1194 

BJNP(1)=CDSIN(Z)/Z  1195 

B0NP(2)=(BJNP(1)-CDC0S(Z))/Z  1196 

ZTI=DCMPLX(-DIMAG(Z),DREAL(Z))  1197 

T1=CDEXP(ZTI)/Z  1198 

T1=DCMPLX (DIMAG(T1 ) ,-DREAl (T1 ) ) 1199 

BHNP(1)=T1  1200 

BHNP(2)=DCMPLX(DIMAG(T1),-DRCAL(T1))*(UD0-1.D0/ZTI)  1201 

ZS0=Z*Z  1202 

TDZ=2.D0/Z  1203 

X=l.n0/ST0PR  1204 

DO  15  N=3,NBF  1205 

XNPH=DFL0AT(N)-.5D0  1206 

XNU=-{XNPH+1.D0)*TDZ  1207 

A1=XNPH*TDZ  1208 

DEN=XNU+1.D0/A1  1209 

F=XNU/(DEN*A1)  1210 

CF=-TDZ  1211 

DO  5 1=2,200  1212 

CF=-CF  1213 

A1=CF*{XNPH+I)  1214 

XNU=A1+1.D0/XNU  1215 

DEN=A1+1. DO/DEN  1216 

F1=XNU/DEN  1217 

F=F*F1  1218 

IF  (DABS(CDABS(F1)-1.D0).LT.1.D-14)  GO  TO  10  1219 

5 CONTINUE  1220 

10  RJNP(N)=F*BJNP(N-1)  1221 

Q=1.D0/(ZSQ*BJNP(N-1))  1222 

BHNP (N )=F*BHNP (N-1 )-DCMPLX {-DIMAG{0 ) ,DREAL (Q ) ) 1223 

IF  (CDABS(BJNP{N)).LT.X.OR.CDABS(BHNP(N)).GT.STOPR)  GO  TO  20  1224 

15  CONTINUE  1225 

N=N-1  1226 

20  IF  {N.LT.5)  PRINT  25,N,Z  1227 

25  FORMAT  {25X, 'ONLY'.IS, ' BESSEL  FUNCTIONS  FOR  Z =MP2D12.4)  1223 

RETURN  1229 

END  1230 

1231 

1232 

SUBROUTINE  SRBF  (A, Y, AD, YD)  1233 

GET  J,  O',  Y AND  Y'  FOR  NEWTON'S  COOLING  FUNCTION  AND  RETURN  1234 
THE  APPROPRIATE  PART  OF  COMPLEX  VALUES  ADJUSTED  FOR  REAL  1235 

VALUE  CALCULATIONS  1236 

IMPLICIT  REALMS  (A-H,0-Z)  1237 

COMMON  /A/NORG,NREG,NRT,NSBF,NMIN,NC,ICODE  1238 

COMMON  /B/FACT(6,25,18),AJ(6,25,18),BY(6,25,18),XLAMDA(25,18),SBDP  1239 
1(6),RH0P(6),CP(6),BP(6),TCP(6),H  1240 

COMMON  /C/AJ1,S,F,R,I  1241 

C0MPLEX*16  BJ,YF,BJD,BYD  1242 

COMMON  /BES/BJ,YF,RJD,BYD,RJ,RY,RJD,RYD  1243 

TP  ^ 1 1^44 

5 CALL  CSBFD(DCMPLX(O.DO,R*F))  1245 

FOR  S<0.  1246 

IF  (NSBF.EQ.2*(NSBF/2))  GO  TO  10  1247 

FOR  S<0.  AND  EVEN  ORDER  BESSEL  FUNCTIONS  1248 

A=DREAL(BJ)  1249 

Y=DIMAG(YF)  1250 
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IF  (ICODE.EO.l)  GO  TO  25  1251 

C=TCP(I)*F  1252 

AD=-C*DIMAG(BJO)  1253 

YO=C*OREAL{ByD)  1254 

GO  TO  25  1255 

FOR  S<0.  AND  ODD  ORDER  BESSEL  FUNCTIONS  1256 

10  A=DIMAG(8J)  1257 

Y=DREALfYF)  1258 

IF  (ICODE.EQ.l)  GO  TO  25  1259 

C=TCP{I)*F  1260 

AD=C*DREAL(BJD)  1261 

YD=-C*niMA6{BYD)  1262 

GO  TO  25  1263 

FOR  S=0.  1264 

15  A=R**(NS8F-1)  1265 

Y=1.D0/R**(NSBF)  1266 

IF  (ICODE.EQ.l)  60  TO  25  1267 

AD=TCP ( I )*(NS8F.l ) *R**(NS8F-2)  1268 

YD-—TCP  ( I )*  (NS8F  )/R**(NSBF+l ) 1269 

GO  TO  25  1270 

FOR  S>0.  1271 

20  CALL  SBFAD(R*F)  1272 

A=RJ  1273 

Y=RY  1274 

IF  (ICODE.EQ.l)  GO  TO  25  1275 

C=TCP(I)*F  1276 

AD=C*RJD  1277 

YD»C*RYO  1278 

25  RETURN  1279 

END  1280 

1281 

1282 

SUBROUTINE  SBFAD(Z)  1283 

SPHERICAL  BESSEL  FUNCTIONS  OF  THE  FIRST  1284 

AND  SECOND  KINDS  AND  THEIR  FIRST  DERIVATIVES  1285 

FOR  REAL  ARGUMENT  1286 

IMPLICIT  REAL*8  (A-H,0-Z)  1287 

COMMON  /A/N0R6,NREG,NRr.NSBF,NMIN,NC,IC0DE  1288 

C0MPLEX*16  BJ.YF.BJD.BYD  1289 

COMMON  /BES/80.YF,BJD,BYD,RJ,RY,RJD.RYD  1290 

COMMON  /C/AJ1,S,F,R,I  1291 

SINZ  = DSIN(Z)/Z  1292 

COSZ  = DC0S(Z)/2  1293 

Y1  = -COSZ  1294 

RY  = Yl/Z  - SINZ  1295 

IF(NSBF.GE»3)  60  TO  12  1296 

iF(NSBF.GT.l)  GO  TO  25  1297 

RJ  = SINZ  1298 

RY=.COSZ  1299 

IF(ICOnE.EQ.l)  GO  TO  55  1300 

RvJD=  COSZ  - SIN7./Z  1301 

RYD  = SINZ  + COSZ/Z  1302 

GO  TO  55  1303 

12  IF  (I.EO.l)  GO  TO  25  1304 

DO  15  M = 3.NSBF  1305 

Y0=Y1  1306 

Y1=RY  1307 
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IF  (DABS(Yl).GT. 1.034)  PRINT  500,NRT,M,NSRF,Z,Y1  1308 


500  FORMAT  ('  SBFAD:  ROOT '.1 3,'  FOR 
1*  Y =’,ni2.4) 

15  RY  - (2*M-3)*Y1/Z  - YO 
25  C - DABS{Z) 

IF(C.GE.3.D0)  GO  TO  30 
RJ  = BES1(NSBF.1,Z) 

GO  TO  35 

30  R.J  = SBFJ(NSBF-1,Z) 

35  IF(ICODE.EQ.l)  GO  TO  55 
IF{NSBF.GT.2)  GO  TO  40 
BJ1=  SINZ 
GO  TO  50 

40  IF(C.GE.3.D0)  GO  TO  45 
BJl  = BESl  (NSBF-2,Z) 

GO  TO  50 

45  BJl  = SBFJ  (NSBF-2,Z) 

50  RJD  = BJl  - NSBF*RJ/Z 
RYD  = Y1  - NSBF*RY/Z 
F/5  RETURN 
END 


FUNCTION  BES1(N,Z) 

IMPLICIT  REAL *8  (A-H,0-Z) 
BES1=DSIN(Z)/Z 
IF  (N.EO.O)  GO  TO  15 
TDZ=2.D0/Z 
11=0 

no  10  M=1,N 

XNUPH=DFL0AT(M)+.5D0 

AO=XNUPH*TDZ 

A1=-(XNUPH+1.D0)*TDZ 

RNUM=A1+1. DO/AO 

RDEN=A1 

COLD=AO*RNUM/RDEN 
CFAC=-TDZ 
DO  5 1=2,200 
CFAC=-CFAC 
A=CFAC*(XNUPH+I) 

RNUM=A+1.D0/RNUM 

RDEN=A+1.D0/RDEN 

C=RNUM/RDEN 

C0LD=C0LD*C 

IF  (DABS(DABS(C)-l.D0).LT.l.D-8) 
5 CONTINUE 
10  BES1=8ES1/C0LD 
15  RETURN 
END 


FUNCTION  S3FJ(N,Z) 
IMPLICIT  REAL *8  (A-H.O-Z) 
Q=O.DO 
P=1.D0 

IF  (N.EO.O)  GO  TO  10 
XN1=N+1 


BFM3,'  OF', 13,'  Z =’,1PD12.4,  1309 

1310 

1311 

1312 

1313 

1314 

1315 

1316 

1317 

1318 

1319 

1320 

1321 

1322 

1323 

1324 

1325 

1326 

1327 

1328 

1329 

1330 

1331 

1332 

1333 

1334 

1335 

1336 

1337 

1338 

1339 

1340 

1341 

1342 

1343 

1344 

1345 

1346 

1347 

1348 

1349 

1350 

1351 

60  TO  10  1352 

1353 

1354 

1355 

1356 

1357 

1358 

1359 

1360 

1361 

1362 

1363 

1364 
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1365 

1366 

1367 

1368 

1369 

1370 

1371 

1372 

1373 

1374 

1375 

1376 

1377 

1378 

1379 

1380 

1381 

1382 

1383 

1384 

1385 

1386 

1387 

1388 

1389 

1390 

1391 

1392 

1393 

1394 

1395 

1396 


Y1  = -COSZ 
YF  = Yl/Z  - SINZ 
IF(NSBF.GE.3)  60  TO  12 
IF(NSBF.GT.l)  GO  TO  25 
BJ  = SINZ 
YF=-COSZ 

IF(ICODE.EQ.l)  GO  TO  55 
BJ0=  COSZ  - SINZ/Z 
BYD  = SINZ  + COSZ/Z 
GO  TO  55 

12  IF  (I.EQ.l)  GO  TO  25 
DO  15  M = 3,NSBF 
Y0=Y1 
Y1=YF 

15  YF  = (2*M-3)*Y1/Z  - YO 

25  IF'C.GE.15.00)  GO  TO  30 
BJ  = BESlC(NSBF-l.Z) 

GO  TO  35 

30  BJ  = SBFJC(NSBF-l.Z) 

35  IF(ICODE.EO.l)  GO  TO  55 
1F(NSBF.GT.2)  60  TO  40 
BJ1=  SINZ 
GO  TO  50 

40  IF(C.6E,15.D0)  GO  TO  45 
BJl  = BES1C(NSBF-2.Z) 


Z2=2.nO*Z 
T=1.D0 

5 T=T*((XN1*XN2)/(F*Z2)) 

Q=Q+T 

IF  {XN2.EQ.1.00)  GO  TO  10 

XN1=XN1+1.D0 

XN2=XN2-1.D0 

F=F+1.D0 

T=-T*f(XNl*XN2)/(F*Z2)) 

P=P+T 

IF  (XN2.E0.1.D0)  GO  TO  10 

XN1=XN1+1.D0 

XN2=XN2-1.D0 

F=F+1.D0 

GO  TO  5 

10  A=Z-N*1.57O7963267948965D0 
SBFJ=(P*OSIN(A)+0*OCOS(A) )/Z 
RETURN 
END 


SUBROUTINE  CSBFD(Z) 

IMPLICIT  C0MPLEX*16  (A-H,0-Z) 

COMMON  /BES/BJ.YF,BJD.BYD 

COMMON  /A/NORG , NREG , NRT , NSBF , NMI N . NC , I CODE 

COMMON  /C/AJ1,S,F,R.I 

REAL*8  C,AJ1,S.F.R 

C * CDABS(Z) 

SINZ  = CDSIN(Z)/Z 
COSZ  = COCOS fZl/Z 


GO  TO  50 

45  BJl  = SRFJC(NSBF-2,7.) 
50  BJD  = BJl  - NSBF*BJ/Z 
BYD  = Y1  - NSBF*YF/Z 
55  RETURN 
END 


FUNCTION  BES1C(N,Z) 

IMPLICIT  C0MPLEX*16  (A-H,0-Z) 

BESIC  = CDSIN(Z)/Z 
IF(N.EQ.O)  GO  TO  15 
BES1C=(BES1C-CDC0S(Z))/Z 
IF  (N.EQ.l)  GO  TO  15 
TDZ  = 2.D0/Z 
DO  10  M = 2,N 

CM  = DCMPLX(DFLOAT(M),O.DO) 

XNUPH  = CM  + .500 
AO  = XNUPHnOZ 
A1  = -(XNUPH  + l.DO)*TDZ 
RNUM  = A1  + 1. DO/AO 
RDEN  = A1 

COLD  = AO*RNUM/RDEN 
CFAC  = -TDZ 
DO  5 I = 2,200 
Cl  = DCMPLX(DFLOAT(I).O.DO) 

CFAC  = -CFAC 
A = CFAC*(XNUPH  + Cl) 

RNUM  = A + l.DO/RNUM 
RDEN  = A + l.DO/RDEN 
C = RNUM/RDEN 
COLD  = COLD*C 

IF(DABS{CDABS(C)-l.D0).LT.l.D-8)  GO  TO  10 
5 CONTINUE 

10  BESIC  = BESIC/COLD 
15  RETURN 
END 


FUNCTION  SBFJC(N,Z) 

IMPLICIT  COMPLEXES  (A-H,0-Z) 
REAL*8  XN1.XN2,F.DREAL, DIMAG 
Q=O.DO 
P=1.D0 

IF  (N.EO.O)  GO  TO  10 

XN1=N+1 

XN2=N 

F=1.D0 

Z2=2.D0*Z 

T=1.D0 

5 T=T*({XN1*XN2)/(F*Z2)) 

0=0+T 

IF  (XN2.EQ.1.D0)  GO  TO  10 

XN1=XN1+1.D0 

XN2=XN2-1.D0 

F=F+1.D0 

T=-T*{(XN1*XN2)/(F*Z2)) 


P=P+T  1479 

IF  (XN2.E0.1.D0)  60  TO  10  1480 

XN1=XN1+1.D0  1481 

XN2»XN2-1.D0  1482 

F=F+1.D0  1483 

60  TO  5 1484 

10  A = 7 - DCMPLX{DFL0AT(N)*1. 570796326794896500,0.00)  1485 

T = (P*CDSIN(A)+Q*COCOS(A))/Z  1486 

IF  (OREAL(Z).EQ.O.OO)  60  TO  17  1487 

SBFJC=T  1488 

60  TO  20  1489 

17  IF  {N.NE.2*(N/2))  60  TO  15  1490 

SBFJC=DCMPLX(DREAL(T),O.DO)  1491 

60  TO  20  1492 

15  SBFJC=DCMPLX(0.D0,DIMA6(T))  1493 

20  RETURN  1494 

ENO  1495 

1496 

1497 

SUBROUTINE  EPROP(F  .ITIS,EPS.SI6)  1498 

IMPLICIT  REAL*8  (A-H,0-Z)  1499 

INTERPOLATE  EPS  ANO  SIGMA  FROM  TABLES  1500 

F FREQUENCY  IN  ME6AHERTZ  1501 

ITIS  TISSUE  TYPE  1502 

1 OENOTES  CEREBROSPINAL  FLUID  1503 

2 DENOTES  BLOOD  1504 

3 DENOTES  MUSCLE  1505 

4 OENOTES  SKIN  OR  DURA  1506 

5 DENOTES  BRAIN  1507 

6 DENOTES  FAT  OR  BONE  1508 

7 DENOTES  YELLOW  BONE  MARROW  1509 

EPS  REAL  PART  OF  DIELECTRIC  CONSTANT  1510 

SIG  CONDUCTIVITY  1511 

DIMENSION  FR(32),EA(32,7),SA(32,7),SA1(128),SA5(96),EA1{128),EA5(9  1512 
16)  1513 

EQUIVALENCE  (SA1,SA),(SA5,SA(1.5)),{EA1,EA),(EA5,EA(1,5))  1514 


DATA  FR/0. 108, . 1259D8, . 1585D8, . 1995D8, .2512D8, . 3162D8, . 3981D8, . 501  1515 
12D8,. 63108,. 794308,. 109,. 1259D9,.1585D9,.1995D9,.2512D9,.3162D9,. 3 1516 
1981D9,.5012D9,.631D9,.7943D9,.1D10, .1259010,. 1585010,. 1995D10,. 251  1517 
12D10,.3162D10,.3981D10,.5012D10,.631D10,.7943D10,.8913D10,.1D11/  1518 

DATA  SAl / . 75D0 , . 762D0 , . 78D0 , . 79800 , . 816D0 , . 8400 , . 876D0 , . 9D0 , . 96D0 ,1519 
1 . 102D1 ,.11401,. 122401 , . 130801 , . 1392D1 , . 1452D1 , . 1524D1 , . 1572D1 , . 160  1520 
18D1 , . 1644D1 , . 174D1 , . 1812D1 , . 1932D1 , . 2064D1 , . 2292D1 ,.261601,. 3084D1  1521 

1.. 3744D1,. 471601,. 64201,. 918D1,. 107602,. 1236D2,.6875D0,.6985D0,. 71  1522 
15D0,.7315D0,. 74800,. 7700,. 803D0,.825D0,.88D0,.935D0,.1045D1,.1122D  1523 

1 1 . .  1 19901 . . 127601 .. 133101 . . 139701 ..144101.. 147401 . .150701 . .159501 , 1524 

1.166101 . . 177101 . . 189201 . .210101 . .239801 . . 282701 . .343201 . .432301 . . 5 1525 
1885D1,. 841501,. 986701,. 113302, .62500,. 63500,. 6500,. 66500,. 68D0,.7D  1526 

10..  7300.. 7500.. 800.. 8500.. 9500.. 10201. .10901.. 11601.. 12101.. 12701,  1527 

1.13101..  13401.. 13701.. 14501. .15101. .16101.. 17201.. 19101.. 21801.. 25  1528 

1701. . 31201.. 39301.. 53501.. 76501.. 89701. .10302.. 531300.. 539300. 0.55  1529 

12500..  565300. 0.57800.. 59500.. 620500.. 637500.. 6800.. 722500.. 807500,  1530 

1.86700..  926500. .98600. .102901. .10801. .111401. .113901. .116501. .1233  1531 

101 . . 128401 . . 136901 . . 146201 . . 162401 . . 185301 . .218501 . .265201 . . 334101  1532 

1 . . 454801 . .650301 . .762501 ..875501/  1533 
DATA  SA5/. 411600,. 418100, 0.428000,. 437900,. 447800,. 46100,. 480700,.  1534 

1493900 . .  526800 . . 559700 . .625600 . .671 700 . . 7 1 7800 . . 763900 . . 796800 . . 83  1535 
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163D0 , .8626D0 , .P324D0 , . 9021D0 , .954800, . 994300 , . 10601 , . 1 1 3301 , . 12580  1536 

11 . .  143601 . . 169201 . . 205501 . .258801 . .352301 . . 503801 . . 590701 . .678301 .1537 
1.220-1,. 2280-1,. 2350-1,. 250-1, .260-1,. 280-1,. 320-1,. 3480-1,. 380-1,  1538 
1.40-1,. <..50-1,. 470-1,. 520-1,. 570-1,. 6280-1,. 690-1,. 740-1,. 810-1,. 8 1539 
180-1,. 960-1,. 10300,. 11300,. 12400,. 13800,. 15400,. 17600,. 20100,. 2360  1540 

10..  27400.. 34200.. 38400.. 436500.. 198D-1,. 20520-1,. 21150-1,. 2250-1,.  1541 
12340-1,. 2520-1,. 2880-1,. 31320-1,. 3420-1,. 360-1,. 38250-1,. 4230-1,. 4 1542 
1680-1 , .5130-1 , . 56520-1 , .6210-1 , .6660-1 , . 7290-1 , . 7920-1 , .8640-1 , .92  1543 
1 70-1 , . 101 700 , . 1 1 1600 , . 124200 , . 138600 , .158400 , . 180900 , . 212400 , . 2466  1544 

100. . 307800. .345600. .392900/  1545 
OATA  EAl/. 19203,. 178203,. 16503,. 151703,. 139303,. 127703,. 117103,. 10  1546 

17203..  987602.. 914402.. 841202.. 771602.. 711602.. 674402.. 658802.. 6480  1547 

12. . 639602 . .632402 . . 62402 . .615602 . .607202 . . 597602 . . 586802 . . 57602 . . 5 1548 

165202..  553202.. 541202.. 526802.. 513602.. 49802.. 489602.. 478802.. 1920  1549 

13..  178203.. 16503.. 151703.. 139303. .127703.. 117103.. 107203.. 987602..  1550 

1914402..  841202.. 771602.. 711602.. 674402.. 658802.. 64802.. 639602.. 632  1551 

1402..  62402.. 615602.. 607202. .597602.. 586802.. 57602.. 565202.. 553202,  1552 

1.541202..  526802.. 513602.. 49802.. 489602.. 478802.. 1603.. 148503.. 1375  1553 

103..  126403.. 116103.. 106403.. 97602.. 89302.. 82302.. 76202.. 70102.. 643  1554 

102..  59302.. 56202.. 54902.. 5402.. 53302.. 52702.. 5202.. 51302.. 50602.. 4 1555 

19802..  48902.. 4802.. 47102.. 46102.. 45102.. 43902.. 42802.. 41502.. 40802  1556 

1..  39902.. 142403.. 132203.. 122403.. 112503.. 103303.. 94702.. 868602.. 79  1557 

14802 . . 732502 . .678202 . .623902 . . 5'23D2 , . 527802 , . 500202 , . 488602 , . 4806  1558 

102. . 474402. .46902. .462802. .456602. .450302. .443202.0.435202. .427202  1559 

1..  419202.. 410302.. 401402.. 390702. .380902.. 369402.. 363102.. 355102/  1560 
OATA  EA5/. 105403,. 977902,. 905402,. 832302, 0.764502,. 700602,. 642702,  1561 

1.58802..  541902.. 501802. 0.461602.. 423402.. 390502.. 370102.. 361502.. 3 1562 

155602..  35102.. 347002. .342402.. 337802. .333202.. 327902.. 32202.. 31610  1563 

12..  310202. .303602.. 29702. .289102. .281802. .273302.. 268702.. 262702..  1564 

13602. . 31802. .27902. .24302. 0.20802.. 17802. . 14802. . 12302. .10'’02, .860  1565 

11..  74501.. 6801.. 6301.. 601.. 5801. .5701. .56501.. 56301. .5620]  56101,  1566 

1.5601..  55901.. 55701.. 55601.. 55401. .55201. .5501.. 54801.. 52Di,. 4901,  1567 

1.47301..  4501.. 32402.. 286202.. 251102. .218702. .187202.. 1602^2,. 13320  1568 

12..  110702.. 91801.. 77401.. 670501.. 61201.. 56701.. 5401.. 52201.. 51301,  1569 

1 . 508501 . .  506701 . . 505801 . . 504901 . . 50401 .. 503101 .. 501301 . . 500401 . . 49  1570 


18601 , .496801 , .49501 , .493201 , .46801 , .44101 , .425701 , .40501/  1571 

FRE0=F  *1.06  1572 

NOX=ITIS  1573 

IF  (FREQ. GE.FR(l). AND. FRE0.LE.FR(32))  GO  TO  10  1574 

PRINT  5 1575 

5 FORMAT  ('0****  FREQUENCY  LIES  OUTSIDE  THE  RANGE  10  - 10000  MHZ  ***  1576 

1*')  1577 

STOP  1578 

10  DO  15  IJ=1,31  1579 

IF  (FRE0.LE.fr ( I J+1))  GO  TO  20  1580 

15  CONTINUE  1581 

20  X=(FRE0-FR(IJ))/(FR(IJ+1)-FR(IJ))  1582 

EPS=EA(IJ,NDX)+(EA(IJ+1,NDX)-EA(IJ,N0X))*X  1583 

SIG=SA(IJ,NDX)+(SA(IJ^I,NDX)-SA(IJ,NDX))*X  1584 

RETURN  1585 

END  1586 

1587 

1588 

SUBROUTINE  PLTCV1(X,Y,XLEN,YLEN,XTL,YTL,NXTL,NYTL,NP,ICRCT,ISYM,  1589 

1 IMM,XMIN,XMAX,YMIN,YMAX,INPLT,LINTYP,SOCH,DELX,DELY,  1590 

2 NDEC)  1591 

WE  ARE  PLOTTING  Y AS  A FUNCTION  OF  X 1592 
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****  THIS  IS  A VARIATION  OF  PLTCRV  TO  PERMIT  SPECIFYING  THE  BLIP  1593 

INTERVAL  AND  THE  NUMBER  OF  DECIMAL  PLACES  AND  CHARACTER  SIZE  FOR  1594 
SCALE  NUMBERS  AND  LABELS.  1595 

X ARRAY  TO  PLOT  ON  X (HORIZONTAL)  AXIS  - DIMENSION  (NP+2)  1596 

Y ARRAY  TO  PLOT  ON  Y (VERTICAL)  AXIS  - DIMENSION  (NP+2)  1597 

XLEN  LENGTH  IN  INCHES  OF  X AXIS  1598 

YLEN  LENGTH  IN  INCHES  OF  Y AXIS  1599 

XTTL  ARRAY  CONTAINING  X AXIS  TITLE  1600 

YTTL  ARRAY  CONTAINING  Y AXIS  TITLE  1601 

NX  NUMBER  OF  CHARACTERS  IN  XTTL  1602 

NY  NUMBER  OF  CHARACTERS  IN  YTTL  1603 

NP  NUMBER  OF  POINTS  TO  PLOT  IN  ARRAYS  X AND  Y 1604 

ICRCT  0 - PLOT  AXES  AND  LINE  PLOT  1605 

1 - PLO'i  LINE  ON  EXISTING  AXES  1606 

ISYM  CODE  (0-13)  TO  SELECT  SYMBOL  TO  MARK  PLOTTED  POINTS  1607 

IMM  0 - GET  SCALE  END  VALUES  BY  SCANNING  X AND  Y ARRAYS  1608 

1 - GET  SCALE  END  VALUES  FROM  INPUT  ARGUMENTS  1609 

XMIN  MINIMUM  VALUE  ON  X AXIS  1610 

XMAX  MAXIMUM  VALUE  ON  X AXIS  1611 

YMIN  MINIMUM  VALUE  ON  Y AXIS  1612 

YMAX  MAXIMUM  VALUE  ON  Y AXIS  1613 

INPLT  0 - DRAW  SCALES  AND  LINE  1614 

1 - GET  MAXIMA  AND  MINIMA  OF  X AND  Y ARRAYS,  NO  PLOT  1615 


LINTYP  MAGNITUDE  GIVES  FREQUENCY  OF  SYMBOLS  - EVERY  LINTYP  PTS.  1616 


=0  - LINE  PLOT.  NO  SYMBOLS  1617 

>0  - LINE  PLOT  WITH  SYMBOLS  1618 

<0  - NO  LINE,  SYMBOLS  ONLY  1619 

SOCH  CHARACTER  HEIGHT  FOR  TITLE  AND  SCALE  (INCHES)  1620 

DELX  FOR  X AXIS,  POSITIVE  VALUE  TO  DEFINE  UNITS  BETWEEN  TIC  1621 

MARKS  (USER  UNITS).  IF  DELX  = 0.,  TIC  MARKS  WILL  BE  1622 
ONE  INCH  APART.  1623 

DELY  FOR  Y AXIS,  POSITIVE  VALUE  TO  DEFINE  UNITS  BETWEEN  TIC  1624 
MARKS  (USER  UNITS).  IF  DELY  = 0.,  TIC  MARKS  WILL  BE  1625 
ONE  INCH  APART.  1626 

NDEC  NUMBER  OF  DECIMAL  PLACES  IN  SCALE  NUMBERS  1627 

>=0  - SPECIFIES  NUMBER  OF  DECIMAL  PLACES  AFTER  1628 

DECIMAL  POINT  1629 

-1  - ROUNDED  INTEGER  DRAWN  1630 

DIMENSION  X(NP),Y(NP),XTL(1),YTL(1)  1631 

IF(ICRCT.EQ.l)  GO  TO  20  1632 

IF(IMM.EQ.l)  GO  TO  10  1633 

XMIN  = 1.E35  1634 

XMAX  = -1.E35  1635 

YMIN  = 1.E35  1636 

YMAX  = -1.E35  1637 

DO  5 I = 1,NP  1638 

XMIN=AMIN1(X(I).XMIN)  1639 

YMIN=AMIN1(Y(I),YMIN)  1640 

XMAX=AMAX1(X(I),XMAX)  1641 

5 YMAX=AMAX1(Y(I),YMAX)  1642 

IF  (INPLT.EO.l)  RETURN  1643 

10  DELVX  = (XMAX-XMIN)/XLEN  1644 

DELVY  = (YMAX-YMIN)/YLEN  1645 

CALL  BAXIS  (0.,0.,XTL,-NXTL. XLEN, 0.,XMIN,0ELVX,DELX, SOCH, NDEC)  1646 

CALL  BAXIS  (0.,0.,YTL,NYTL, YLEN, 90., YMIN, DELVY.DELY, SOCH, NDEC)  1647 

20  IF(ISYM,LT.0.0R.ISYM.GT.13)  ISYM  = 1 1648 

X(NP+1)  = XMIN  1649 


Y(NP+1)  = YMIN  1650 

X{NP+2)  = (XMAX-XMIN)/XLEN  1651 

Y(NP+2)  = (YMAX-YMIN)/YLEN  1652 

CALL  LINE(X,Y,NP,1.LINTYP,ISYM}  1653 

RETURN  1654 

END  1655 

1656 

1657 

SUBROUTINE  BAXIS  (XPAGE.Y?AGE.IBCD,NCHAR,AXLEN, ANGLE, FIRSTV.DELTAV  1658 
1,DELTIC,S0C::,NDEC)  1659 

THIS  SUBROUTINE  IS  AN  EXTENSION  OF  THE  CALCOMP  'AXIS'  ROUTINE  1660 
TO  ALLOW  THE  USER  TO  SPECIFY  THE  SIZE  OF  CHARACTERS,  THE  1661 

DISTANCE  BETWEEN  TIC  MARKS  AND  THE  NUMBER  OF  DECIMAL  PLACES  IN  1662 
THE  SCALE  NUMBERS.  1663 

XPAGE  - X COORDINATE  OF  AXIS  STARTING  POINT  (INCHES)  1664 

YPAGE  - Y COORDINATE  OF  AXIS  STARTING  POINT  (INCHES)  1665 

IBCD  - ARRAY  WITH  AXIS  TITLE  1666 

NCHAR  - NUMBER  OF  CHARACTERS  IN  AXIS  TITLE  1667 

<0  - ALL  NOTATION  ON  CLOCKWISE  SIDE  OF  AXIS  1668 

>0  - ALL  NOTATION  ON  COUNTERCLOCKWISE  SIDE  1669 

AXLEN  - AXIS  LENGTH  (INCHES)  (MUST  BE  POSITIVE)  1670 

ANGLE  - ANGLE  (POSITIVE  OR  NEGATIVE)  AT  WHICH  AXIS  IS  DRAWN  1671 
(DEGREES)  1672 

FIRSTV  - STARTING  VALUE  (MAX  OR  MIN)  OF  AXIS  AT  FIRST  TIC  1673 
(USER  UNITS)  1674 

DELTAV  - INCREMENT  OR  DECREMENT  VALUE  ASSOCIATED  WITH  ONE  1675 
INCH  ON  AXIS  (USER  UNITS)  1676 

DELTIC  - POSITIVE  VALUE  TO  DEFINE  UNITS  BETWEEN  TIC  MARKS  1677 
(USER  UNITS)  IF  DELTIC  = 0.,  TIC  MARKS  WILL  BE  1678 
ONE  INCH  APART.  1679 

SOCH  - CHARACTER  HEIGHT  FOR  TITLE  AND  SCALE  (INCHES)  1680 

NDEC  - NUMBER  OF  DECIMAL  PLACES  IN  SCALE  NUMBERS  1681 

>=0  - SPECIFIES  NUMBER  OF  DECIMAL  PLACES  AFTER  1682 
DECIMAL  POINT  1683 

-1  - ROUNDED  INTEGER  DRAWN  1684 

DIMENSION  IBCD(l)  1685 

IF  (AXLEN.GT.0..AND.DELTIC.GE.0..AND.NDEC.LE.9)  GO  TO  10  1686 

PRINT  5,AXLEN,DELTIC,NDEC  1687 

5 FORMAT  ('0****  BAXIS  ERROR:  AXLEN  =',1PD15.7,'  DELTIC  =',D15.7,'  1688 
1 NDEC  =',I5,'  ****')  1689 

STOP  1690 

10  IF  (NDEC.LT.-1)N0EC=-1  1691 

AIR=3.1415927*ANGLE/180.  1692 

CA=COS(AIR)  1693 

SA=SIN(AIR)  1694 

DRAW  AXIS  LINE  1695 

CALL  PLOT (XPAGE, YPAGE, 3)  1696 

CALL  PL0T(XPAGE+AXLEN*CA,YPAGE+AXLEN*SA,2)  1697 

FSTV^FIRSTV  1698 

DEI.V=DELTAV  1699 

A=AMAX1(ABS(FSTV),ABS(FSTV+DELV*AXLEN))  1700 

M=AL0G10(A)  1701 

IF  (A.LT..1)M=M-1  1702 

TM=10.*^  1703 

DTIC=ABS(DELTIC/DELV)  1704 

IF  (DELTIC. E0.O)DTIC=l.O  1705 

DELV=DELV/TM  1706 
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FSTV=FSTV/TM  1707 

Xl=S0CH/2.  1708 

TICH=X1  1709 

IF  (NCHAR.LT.O)TICH=-TICH  1710 

XT=-TICH*SA  1711 

YT=TICH*CA  1712 

COMPUTE  POSITION  OF  AXIS  SCALE  NUMBERS  RELATIVE  TO  TIC  1713 

MARKS  AND  ADJUST  FOR  NUMBER  OF  DECIMAL  POINTS  1714 

FN=X1  1715 

IF  (NDEC.GE.0)FN=FN*(2+NDEC)  1716 

FN=FN-.429*X1  1717 

XN=1.4*XT-FN*CA  1718 

YN=1.4*YT-FN*SA  1719 

IF  (NCHAR.GT.O)  GO  TO  20  1720 

FOR  TICS  ON  CLOCKWISE  SIDE  OF  AXIS,  NUMBERS  MUST  BE  MOVED  1721 

AWAY  FROM  AXIS  BY  ONE  CHARACTER  WIDTH  1722 

XN=XN+2.*XT  1723 

YN=YN+2.*YT  1724 

20  XTIC=XPAGE  1725 

YTIC=YPAGE  -"'-S 

DX=0TIC*CA  i72V 

DY=DTIC*SA  -728 

FPN=FSTV  729 

DTIC=DTIC*DELV  1730 

X=.571*S0CH-FN  1731 

IL=0  1732 

LOOP  TO  DRAW  TICS  AND  SCALE  NUMBERS  1733 

25  CALL  PL0T(XTIC,YTIC.3)  1734 

CALL  PL0T(XTIC+XT,YTIC+YT,2)  1735 

IF  (IL.EQ.O.AND.NDEC.GE.O)  60  TO  30  1736 

X=0.  1737 

IF  (FPN.LT.0.)X=X1  1738 

30  CALL  NUMBER(XTIC+XN-X*CA,YTIC+YN-X*SA,SOCH,FPN,ANGLE,NDEC)  1739 

XTIC=XTIC+DX  1740 

YTIC=YTIC>DY  1741 

FPN=FPN+DTIC  1742 

ALEN=(XTIC-XPAGE-DX*.5)/CA  1743 

IF  (ALEN.GT.AXLEN)  GO  TO  45  1744 

IL=IL+1  1745 

IF  (IL.LE.lOO)  GO  TO  25  1746 

PRINT  40  1747 

40  FORMAT  ('0****  8AXIS  ERROR:  MORE  THAN  100  TIC  MARKS  ****')  1748 

STOP  1749 

CENTER  AXIS  TITLE  AND  PLOT  IT  1750 

45  IL=!ABS(NCHAR)  1751 

IF  (M.NE.0)IL=IL+4  1752 

X=rL*SOCH  1753 

HTL=(AXLEN-X)/2.  1754 

XN=XPAGE+4.6*XT+HTL*CA  1755 

YN=VPAGE+4.6*YT+HTL*SA  1756 

IF  (NCHAR.GT.O)  GO  TO  50  1757 

LEAVE  ROOM  FOR  TITLE  CHARACTERS  ON  CLOCKWISE  SIDE  OF  AXIS  1758 
XN=XN+2.*XT  1759 

YN=YN+2.*YT  1760 

50  CALL  SYMBOL(XN,YN,SOCH,IBCD,ANGLE,IABS(NCHAR))  1761 

IF  (M.EQ.O)  GO  TO  55  17t,2 

ADD  SCALE  FACTOR  1763 
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CALL  SYMB0L(999.,999.,S0CH,'  *10' .ANGLE, 4)  1764 

XN=XN+X*CA-X1*SA  1765 

YN=YN+X*SA+1.5*X1*CA  1766 

CALL  NUMBER(XN,YN,X1,FL0AT(M).ANGLE,-1)  1767 

55  RETURN  1768 

END  1769 

1770 

1771 

SUBROUTINE  CNTRP1(X,NR0W,Y,NC0L,D,NLEV1,NSYM1.IFL)  1772 

DIMENSION  X{NR0W),Y(NC0L),D(NR0W.NC0L),FLEV(10),XST{60),YST(60)  1773 

INTEGER*2  IFL(NR0W,NC0L),IST(6O),JST(60)  1774 

NLEV=NLEV1  1775 

IF  (NLEV.LT.1)NLEV=1  1776 

IF  {NLEV.GTcl0)NLEV=10  1777 

NSYM=NSYM1  1778 

IF  (NSYM.LE.O)NSYM=NROW*NCOL  1779 

AXLEN=6.  1780 

SCALE  THE  DATA  FOR  THE  COORDINATE  AXES  1781 

7MAX=-1.E38  1782 

ZMIN=1.E38  1783 

XMAX=-1.E38  1784 

XMIN=1.E38  1785 

DO  5 I=1,NR0W  1786 

IF  (X(I).GT.XMAX)XMAX=X(I)  1787 

IF  (X(I).LT.XMIN)XMIN=X(I)  1788 

DO  5 J=1,NC0L  1789 

IF  (D(I,J).GT.ZMAX)ZMAX=0(I.J)  1790 

5 IF  (0(I,J).GT.O..AND.D(I,J).LT.ZMIN)ZHIN=D(I,J)  1791 

YMAX=-1.E38  1792 

YMIN^l.E+38  1793 

DO  10  J=1,NC0L  1794 

IF  (Y{J),6T.YMAX)YMAa=Y(J)  1795 

10  IF  (Y(J).LT.YMIN)YMIN*Y(J)  1796 

PRINT  15,XMIN,XMAX,YMIN,YMAX,ZMIN,ZMAX  1797 

15  FORMAT  ('OX  RANGE' ,1P2E12.4, ' Y RANGE' ,2E12.4, ' Z RANGE', 2 1798 
1E12.4)  1799 

XFAC=AXLEN/(XMAX-XMIN)  1800 

YFAC=AXLEN/(YMAX-YMIN)  1801 

CDIF1=(ZMAX-ZMIN)/(2*NLEV)  1802 

IL=-AL0G10(CDIF1)+1.  1803 

20  T=10.**IL  1804 

ICDIF=5*( (IFIX{CDIFl*T)+2)/5)  1805 

S=(ZMIN+ZMAX-FL0AT{2*NLEV*ICDIF)/T)/2.  1806 

IS=0  1807 

IF  (S.NE.0.)IS=5*(IFIX(S*T+2.5*S/ABS{S))/5)  1808 

T1=FL0AT(IS+ICDIF)/T  1809 

CDIF=FL0AT(2*ICDIF)/T  1810 

S=T1+CDIF*{NLEV-1)  1811 

Sl=CDIF*.l  1812 

IF  (ZMIN.LT.Tl-SlcAND.ZMAX.GT.S+Sl.AND.ZMAX.LT.S+CDIF)  GO  TO  25  1813 

IL=IL+1  1814 

GO  TO  20  1815 

25  FLEV(1)=T1  1816 

IF  INLEV. EQ.l)  GO  TO  35  1817 

DO  30  K=2.NLEV  1818 

30  FLEV{K)=T1+FL0AT(2*ICDIF*{K-1))/T  1819 

35  AXLPl=AXLEN+.5  1820 
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AXLP2=AXLEN+.75 

RSQ=X(1)**2 

AX2=AXLEN/2. 

AX2S=(.985*AX2)**2 
NR0WM1=NR0W-1 
NC0LM1=NC0L-1 
DO  40  K=1,NLEV 
S=FLEV(K)+.001*CDIF 
T=FLEV(K)-.001*CDIF 
DO  40  I=1,NR0W 
DO  40  J=1,NC0L 

40  IF  (D(I,J).LT.S.AND.D(I,J).GT.T)D(I,J)=S 
DO  380  K=1,NLEV 
F=FLEV(K) 

IEND=0 

DO  150  I=1,NR0WM1 
DO  150  J=1,NC0LM1 
IFL(I,J)=0 
DIJ=D(I,J) 

DI1J=D(I+1,J) 

DIJ1=D(I,J+1) 

DI1J1=D(I+1,J+1) 

IF  (DIJ.GT.F.OR.OIJl.LT.F)  GO  TO  85 

T=DI1J1-F 

A=DIJ 

B=DI1J1 

45  IF  (I.GT.l)  GO  TO  60 
IF  (A.GT.O.)  GO  TO  50 
YC=SQRT(RSQ-X(I)**2) 

IF  (Y(J+1).LT.0.)YC=-YC 
GO  TO  55 

50  S=(F-DIJ)/(DIJ1-DIJ) 
YC=Y(J)+S*(Y(J+1)-Y(J)) 

55  IEND=IEND+1 
YST{IEND)=YC 
XST(IEND)=X(I) 

IST(IEND)=0 

JST(IEND)=J 

60  IF  (T.GT.O.)  GO  TO  80 

65  IF  (J.LT.NCOLHl)  GO  TO  80 
IF  (B.GT.O.)  GO  TO  70 
XC=SQRT(RSQ-Y{J+1)**2) 

IF  (X{I+1).LT.0.)XC=-XC 
GO  TO  75 

70  S=(F-DIJ1)/(DI1J1-DIJ1) 
XC=X(I)+S*(X{I+1)-X(I)) 

75  IEN0=IEND+1 
XST(IENDV-’XC 
YST(IEND)=Y(J+1) 

IST{IEND  =I 
JST(IEND)=NCOL 

80  IFL(I,J)=1 
GO  TO  95 

85  IF  (DIJ.LT.F.OR.DIJl.GT.F)  GO  TO  90 
T=F-DI1J1 


GO  TO  45 
90  B=DIJ1 

IF  (DIJ.LT.F.AND.DIIJI.GT.F)  GO  TO  65 
B=DI1J1 

IF  (DIJ.GT.F.AND.DIIJI.LT.F)  GO  TO  65 
95  IF  (DIJ.GT.F.OR.DIIJ.LT.F)  GO  TO  140 
T=DI1J1-F 
A=DIJ 
B=DI1J1 

100  IF  (J.GT.l)  GO  TO  115 
IF  (A.GT.O.)  GO  TO  105 
XC=SQRT(RSQ-Y{J)**2) 

IF  (X{I+1).LT.0.)XC=-XC 
GO  TO  110 

105  S=(F-DIJ)/(DI1J-DIJ) 

XC=X(I)+S*(X(I+1)-X(I)) 

110  IEND=IEND+1 
XST{IEND)=XC 
YST(IEND)=Y{J) 

IST(IEND)=I 

JST(IEND)=0 

115  IF  (T.GT.O.)  GO  TO  135 
IFL(I,J)=IFL(I,J)+2 
120  IF  (I.LT.NROWMl)  GO  TO  135 
IF  (B.GT.O.)  GO  TO  125 
YC=SQRT(RSQ-X{I+1)**2) 

IF  (Y(J+1).LT.0.)YC=-YC 
GO  TO  130 

125  S=(F-DHJ)/(DI1J1-DI1J) 

YC=Y{J)+S*(Y(J+1).Y(J)) 

130  IEND=IEND+1 
YST(IEND)=YC 
XST{IEND)=X(I+1) 

IST(IEND)=NROW 

JST(IEND)=J 

135  IF  (IFL(I,J).EQ,0)IFL(I,J)=1 
GO  TO  150 

140  IF  (DIJ.LT.F.OR.DIIJ.GT.F)  GO  TO  145 
T=F-DI1J1 
A=DnJ 
B=A 

GO  TO  100 
145  B=DI1J 

IF  (DIJ.LT.F.AND.DIIJI.GT.F)  GO  TO  120 
B=DI1J1 

IF  (DIJ.GT.F.AND.DIIJI.LT.F)  GO  TO  120 
150  CONTINUE 

155  IF  (lEND.EQ.O)  GO  TO  160 

SET  UP  TO  PLOT  NEXT  CONTOUR  FROM  EDGE  OF  GRID 
I=IST(1) 

J=JST(1) 

IOLD=I 

jnLD=J 

CALL  PLOT ( ( XST ( 1 ) -XMI N ) *XFAC , ( YST ( 1 ) ~ YMI N ) *y FAC , 3 ) 
IF  (I. EO. 0)1=1 
IF  (J.E0.0)J=1 
IF  (I,EQ.NR0W)I=NR0WM1 


IF  (J.EQ.NC0L)J=NC0LM1  1935 

ISTC=1  1936 

GO  TO  180  1937 

ALL  CONTOURS  THAT  LEAVE  GRID  HAVE  BEEN  DRAWN  1938 

SET  UP  Tc  PLOT  NEXT  CONTOUR  THAT  DOES  NOT  LEAVE  GRID  1939 

160  11=1  1940 

165  Jl=l  1941 

170  F (IFL(I1,J1).NE.0)  GO  TO  175  1942 

J1=J1+1  1943 

IF  (Jl.LT.NCOL)  60  TO  170  1944 

11=11+1  194b 

IF  (Il.LT.NROW)  GO  TO  165  1946 

GO  TO  375  1947 

175  ISTC=0  1948 

1=11  1949 

J=J1  1950 

180  ISYM=NSYM-1  1951 

FIND  ENDS  OF  LINES  IN  UPPER  LEFT  TRIANGLE  1952 

185  DIJ=D(I,J)  1953 

DI1J=D(I+1,J)  1954 

DIJ1=D(I.J+1)  1955 

0I1J1=D(I+1.J+1)  1956 

IF  (DIJ.GT.F.OR.DIJl.LT.F)  GO  TO  220  1957 

T=DI1J1-F  1958 

A=DIO  1959 

B=DI1J1  1960 

190  IF  (A.GT.O.)  GO  TO  195  1961 

YC=SQRT(RSQ-X(I)**2)  1962 

IF  (Y(J+1).LT.0.)YC=-YC  1963 

GO  TO  200  1964 

195  S={F-DIJ)/(DIJ1-DIJ)  1965 

YC=Y(J)+S*(Y(J+1)-Y(J))  1966 

200  XC=X(I)  1967 

IB=I-1  1968 

JB=.J  1969 

IF  (T.GT.O,)  GO  TO  250  1970 

IF  (IFL(I,J).E0.2)  GO  TO  250  1971 

IF  (B.GT.O)  GO  TO  205  1972 

XC1=S0RT(RSQ>Y{J+1)**2)  1973 

IF  (X{I+1).LT.0.)XC1=-XC1  1974 

GO  TO  210  1975 

205  S=(r-DIJl)/(OIlJl-OIJl)  1976 

XC1=X(I)+S*{X(I+1)-X(I))  1977 

210  YC1=Y(J+1)  1978 

IE=I  1979 

JE=J+1  1980 

IF  (lOLD.EQ.IB.AND.JOLD.EQ.JB)  GO  TO  215  1981 

IF  (lOLD.NC.IE.OR.JOLD.NE.JE)  GO  TO  250  1982 

215  IFL(I,J)=IFL(I,J)-1  1983 

GO  TO  310  1984 

220  IF  (DIJ.LT.F.OR.DIJl.GT.F)  GO  TO  225  1985 

T=F-Dnul  1986 

A=DIJ1  1987 

B=A  1988 

GO  TO  190  1989 

225  IF  (DIJ.6T.F.0R.DI1J1.LT.F)  GO  TO  24b  1990 

IF  (DIJUGT.O.)  GO  TO  235  1991 
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230  XC=SQRT(RS0-Y(J+1)**2) 

IF  {X(1+1).LT.0.)XC=-XC 
GO  TO  240 

235  S={F-DIJ1)/(DI1J1-DIJ1) 

XC=X(I)+S*(X(I+l)-XfI)) 

240  vc=Y(J+l) 

IB=I 
JB=J+1 
GO  TO  250 

245  IF  (DIJ.LT.F.OR.DIIJI.GT.F)  GO  TO  250 
IF  (OIIJI.GT.O.)  GO  TO  235 
GO  TO  230 

FIND  ENDS  OP  LINES  IN  LOWER  RIGHT  TRIANGLE 
250  IF  (DIJ.GT.F.OR.DIIJ.LT.F)  GO  TO  290 


T=DI1J1-F 
A=DIJ 
B=DI1J1 

255  IF  (A.GT.O.)  GO  TO  260 

XC1=SQRT(RSQ-Y(J)**2)  2C 

IF  (X{I+1).LT.0.)XC1=-XC1  2C 

GO  TO  265  2C 

260  S=(F-DIJ)/(DI1J-DIJ)  2C 

XC1=X(I)+S*{X(I+1)-X(I))  2( 

265  IF  (T.GT.O)  GO  TO  285  2C 

IF  (IFL{I.J).LT.2)  GO  TO  310  2C 

XC=XC1  2( 

YC=Y(J)  2( 

IB=I  2( 

2( 

IFL(I,J)=IFL(I,J)-2  2( 

270  IF  (B.GT.O.)  GO  TO  275  2( 

YCi=SQRT(RSQ-X(I+l)**2)  2( 

IF  (Y(J+1).LT.0)YC1=-YC1  21 

GO  TO  280  2( 

275  S=(F-DIlJ)/(OU,)).-DIlw)  2( 

YCl=Y(J)+S*(Y(J+})-Y(J))  2( 

280  XCl=xn+l)  21 

IEuI+1  21 

JE=»1  21 

GO  TO  310  2i 


IE=I 
-JE=J-1 
IFL(I,J)=0 
GO  TO  310 

290  IF  {DIJ.LT.F.OR.DIivl.GT.P)  GO  TO  295 
T=F-DiiJl 
A=0I1J 
B=A 

GO  TO  255 

295  IF  (DIO.GT.F.CP.DJMl.LT.F)  GO  TO  305 
B=0I1J 

300  IFL{i.J)=0 
GO  TO  270 
305  B^nilJl 

IF  (DiJcGT.F.AND.OIlul.LT.F)  GO  TO  300 
310  IF  (ISTC.NE.O)  GO  TO  320 
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?iAn  FIRST  SEGMENT  OF  NEW  CONTOUR 

2049 

CALI  Pl.OT ( ( XC 3.-XMI N )*XFAC , ( YCl-YMI N ) *YFAC  ,3) 

2050 

ISTC=: 

2051 

315  PX=(XC-XMIN)*XFAC 

2052 

PY=(YC-YMIN)*YFAC 

2053 

IOLD=I 

2054 

J0LD=J 

2055 

I=IB 

2056 

J--JB 

2057 

GO  TO  340 

2058 

MATCH  CURRENT  PEN  POSITION  TO  ONE  END  OF 

NEW  LINE  SEGMENT 

2059 

320  IF  (lOLD.EQ.IB.AND.JOLD.EQ.JB)  GO  TO  335 

2060 

IF  (lOLD.EQ.IE.AND.JOLD.EQ.JE)  GO  TO  315 

2061 

PRINT  330,I,J,I0LD,J0LD,IB,0B,IE,JE,0IJ,DIJi,DIlJ,DIlJl 

2062 

330  FORMAT  {'-LOGIC  ERROR:  AT', 213/  FROM', 213,' 

TO', 213,'  OR', 213, 4 

2063 

1F8.4) 

2064 

STOP 

2065 

335  PX=(XC1-XMIN)*XFAC 

2066 

PY=(YC1-YMIN)*YFAC 

2067 

IOLD=I 

2068 

J0LD=0 

2069 

I=IE 

2070 

J=JE 

2071 

PLOT  LINE  SEGMENT 

2072 

340  R1SQ=(PX-AX2)**2+(PY-AX2)**2 

2073 

IF  (R1SQ.GT.AX2S)  GO  TO  345 

2074 

ISYM=ISYM+1 

2075 

IF  (ISYM.LT.NSYM)  GO  TO  345 

2076 

CALL  SYMB0L(PX,PY,.07,K,0.,-2) 

2077 

ISYM=0 

2078 

60  TO  350 

2079 

345  CALL  PL0T(PX,PY,2) 

2080 

DETERMINE  WHETHER  CONTOUR  HAS  ENDED 

2081 

350  IF  (I.EO.O.OR.I.EO.NROW.OR.J.EQ.O.OR.O.EQ.NCOL)  GO  TO  355 

2082 

IF  (IFL(I,J).NE,0)  GO  TO  185 

2083 

IF  (lEND.EQ.O)  GO  TO  170 

2084 

REMOVE  END  POINTS  OF  LAST  CONTOUR  FROM  TABLE 

2085 

355  11=0 

2086 

IF  (lEND.EQ.l)  GO  TO  370 

2087 

DO  365  L=2,IEND 

2088 

IF  (I.EQ.IST(L).AND.J.EQ.JST(L))  60  TO  365 

2089 

11=11+1 

2090 

XST(I1)=XST(L) 

2091 

YST(I1)=YST(L) 

2092 

IST(I1)=IST(L) 

2093 

JST(I1)=JST(L) 

2094 

365  CONTINUE 

2095 

370  IEND=I1 

2096 

GO  TO  155 

2097 

PUT  SYMBOL  AND  LEVEL  ON  PLOT 

2098 

ALL  CONTOURS  AT  THIS  LEVEL  HAVE  BEEN  DRAWN 

2099 

375  FLK=FL0AT(K-1)*.6 

2100 

CALL  SYMB0L(AXLP1,FLK+.08,.14,K,0.,-1) 

2101 

CALL  FNUM(AXLP2,FLK,FLEV(K),2,0.,.14) 

2102 

380  CONTINUE 

2103 

***  PLOT  AXES 

2104 

CALL  PL0T(AX2,AXLEN,3) 

2105 

CALL  PL0T(AX2,0.,2)  2106 

CALL  PL0T(0.,AX2,3)  2107 

CALL  PL0T(AXLEN,AX2,2)  2108 

***  PLOT  CIRCLE  AROUND  CONTOURS  2109 

DTH=2. *3. 1415927/288  2110 

THETA=DTH  2111 

no  385  K=l,288  2112 

R=AX2*(1,+C0S (THETA))  2113 

P=AX2*(l.fSIN(THETA))  2114 

THETA=THETA+DTH  2115 

385  CALL  PLCT(R,P,2)  2116 

RETURN  2117 

END  2118 

2119 

2120 

SUBROUTINE  rNUM(XPAGE,YPAGE»FPN,ND, ANGLE, HEIGHT)  2121 

EDIT  A FLOATING  POINT  NUMBER  FOR  THE  PLOTTER  2122 

XPAGE  X COORDINATE  OF  STARTING  POINT  (INCHES)  2123 

YPAGE  Y COORDINATE  OF  STARTING  POINT  (INCHES)  2124 

FPN  NUMBER  TO  BE  PLOTTED  2125 

ND  IF  10.**-ND  <=  FPN  < 10.**ND,  THE  NUMBER  IS  PLOTTED  2126 
WITHOUT  EXPONENT  2127 

ANGLE  ANGLE  AT  WHICH  NUMBER  IS  PLOTTED  (DEGREES)  2128 

HEIGHT  CHARACTER  SIZE  (INCHES)  2129 

DIMENSION  8(6)  2130 

DATA  B/. 999999,. 99999,. 9999,. 999,. 99,. 9/  2131 

X=ABS(FPN)  2132 

M=ND  2133 

IF  (X.NE.O.)  GO  TO  5 2134 

PLOT  ZERO  2135 

CALL  NUMBER(XPAGE,YPAGE,HEIGHT,X,ANGLE,1)  2136 

RETURN  2137 

5 N=ALOG10(X)  2138 

IF  (X.LT.1.)N=N-1  2139 

X=X*10.**(-N)  2140 

T=X/10.-X*l.E-7  2141 

no  10  J=l,6  2142 

T1=T-INT(T)  2143 

IF  (Tl.LE.O)  GO  TO  15  2144 

IF  (Tl.GE.B(J))  GO  TO  15  2145 

10  T=T*10.  2146 

15  J=J-1  2147 

T=ABS(FPN)  2148 

IF  (T.GE.10.**M)  GO  TO  20  2149 

IF  (T+.5*10.**(N-6).LT.10.**(-M))  GO  TO  20  2150 

PLOT  NUMBERS  WHICH  DO  NOT  NEED  EXPONENTS  2151 

M=J-N-1  . 2152 

IF  (M.LT.1)M=1  2153 

IF  (M.GT.9)M=9  2154 

CALL  NUMBER(XPAGE. YPAGE, HEIGHT, FPN, ANGLE,M)  2155 

RETURN  2156 

PLOT  NUMBERS  WITH  EXPONENTS  2157 

20  IF  (J.GT.l)  GO  TO  25  2158 

X=l.  2159 

N=N+1  2160 

J=2  2161 

25  IF  (FPN.LT.O.)X=-X  2162 
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CALL  NUMBER (XPAGE , YPA6E .HE IGHT, X .ANGLE .0-1) 
CALL  SYMBOL (999. .999. .HE IGHT. 3H*10. ANGLE . 3) 
Xl=HEIGHT/2 

A=3.1415927*ANGLE/180. 

SA=SIN(A) 

CA=COS(A) 

NC=J+4 

IF  (FPN.LT.0.)NC=NC+1 

S=NC*HEIGHT 

PX=XPAGE+S*CA-X1*SA 

PY=YPAGE+S*SA+1.5*X1*CA 

CALL  NUMBER (PX , PY . XI . FLOAT{N ) .ANGLE . -1 ) 
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DEPARTMENT  OF  THE  AIR  FORCE 
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MEMORANDUM  FOR  DTIC-OCQ 

ATTf^t  LARRY  DOWNING 

8725  JOHN  J.  KINGMAN  ROAD,  SUTTE  0944 

FORT  BELVOIR,  VA  22060-6218 

FROM:  AFIOH/DOBP  (STINFO) 

2513  Kennedy  Circle 

Brooks  City-Base  TX  78235-5 1 1 6 

SUBJECT;  Changing  the  Distribution  Statement  on  a Technical  Report 


This  letter  documents  the  requirement  for  DTIC  to  change  the  distribution  statement  from  “B”  to  “A” 
(Approved  for  public  release;  distribution  is  unlimited.)  on  the  following  technical  report:  AD  Number 
ADB071 126,  SAM-TR-82-22,  A Computer  Model  Predicting  the  Thermal  Response  to  Microwave 
Radiation. 


If  additional  information  or  a corrected  cover  page  and  SF  Form  298  arc  required  please  let  me  know.  You 
can  reach  me  at  DSN  240-6019  or  my  e-mail  address  is  sherrv.mathews@brooks.af.mil. 


Thank  you  for  your  assistance  in  making  this  change. 


SHERRY^.  MATHEWS 
AHOH  STINFO  Officer 


